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CHAPTER  1 


INTRODUCTION 


1 . I General 

Elastic  analysis  produces  satisfactory  results  when  the 
loading  develops  stress  below  the  elastic  limit.  But  it  obviously  begins 
to  experience  trouble  in  predicting  the  structural  behavior  once  the 
yield  stress  has  been  exceeded.  For  intersecting  cylinders,  a stress 
concentration  exists  in  the  vicinity  of  the  intersection  curve.  It  is 
therefore  impractical  or  at  least  uneconomical  to  let  this  concentration 
control  the  design  while  retaining  the  same  allowable  stress  throughout 
the  whole  structure  unless  that  concentrated  stress  critically  determines 
the  failure  load.  Actually,  to  tolerate  a small  amount  of  plastic 
deformation  in  the  region  of  high  stress  gradient  helps  the  material  to 
accommodate  the  imposed  distortion  pattern  and  smooths  out  the  stress 
concentration  as  long  as  the  material  is  ductile  and  no  fatigue  crack 
occurs.  For  structures  operating  in  a high  pressure  state  or  for  large 
diameter  intersecting  cylinders,  this  allowance  saves  material  and  serves 
as  a safety  valve.  Design  procedures  based  on  this  principle  have  been 
formalized,  for  example,  the  ASME  Pressure  Vessel  Code  [1]  allows  self- 
equilibrating  thermal  stresses  calculated  by  elastic  procedures  to  be  up 
to  twice  the  value  of  the  yield  stress.  However,  the  accurate  and  detailed 
determination  of  the  stresses  in  the  vicinity  of  the  intersection  region 
of  the  cylinders  would  be  of  little  value  unless  the  designer  recognized 
the  significance  of  those  stresses  in  relation  to  failure.  A better 
understanding  of  the  post-elastic  behavior  and  the  possibility  of  achieving 
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a more  reliable  design  of  the  cylindrical  intersections  can  be  obtained 
from  a plastic  analysis. 

Cyl inder-to-cyl inder  intersections  are  a very  common  occurrence 
in  many  industrial  applications  such  as  boilers,  pressure  vessels,  pipe 
connections,  etc.  However,  until  only  a few  years  ago  most  of  the 
research  investigations  reported  in  the  literature  were  limited  to 
experimental  work.  Recently  analytical  treatment  of  this  subject  area 
has  been  given  some  attention.  Most  of  this  recent  work  is  still 
incomplete  [2]. 

Difficulties  in  obtaining  analytical  evaluations  of  the  stress 
distributions  in  the  disturbed  regions  near  the  intersection  of  comparable 
size  shells  originally  stemmed  from  the  complicated  geometrical  shape  of 
the  intersection  line.  The  intersection  curve  of  the  middle  surfaces  of 
the  cylinders  is  neither  rotational  symmetric  nor  on  a plane  curve  but 
rather  is  a spacial  curve.  Early  efforts  required  one  cylinder  to  be  of 
a much  smaller  diameter  in  comparison  to  the  cylinder  that  it  is  inter- 
secting so  that  the  intersection  curve  could  be  approximated  by  a circle, 
thus  simplifying  the  problem.  Besides,  the  sharp  discontinuities  of 
curvatures  across  the  intersection  curve  function  as  a stress  raiser. 
Therefore,  the  presence  of  the  stress  concentration  is  inevitable  and, 
as  a consequence,  constitutes  a major  consideration  in  the  design. 

With  the  aid  of  high-speed  digital  computers,  numerical 
solutions  are  now  playing  a significant  role  in  obtaining  solutions  to 
engineering  applications.  During  the  past  decade,  the  development  of  the 
finite  element  method  has  increased  markedly  the  capability  of  engineeting 
problem  solving.  Many  complicated  design  problems  which  were  considered 
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unfeasible  to  a realistic  analysis  several  years  ago  can  now  be  solved 
almost  routinely  by  using  the  finite  element  method.  The  method  provides 
a powerful  tool  to  attack  shell  structures  and  has  been  applied  lately  to 
evaluate  the  stress  distribution  for  intersecting  cylinders  [3].  By 
subdividing  the  whole  structure  into  a finite  number  of  regions,  referred 
to  as  the  "elements",  it  has  the  advantage  of  being  able  to  adjust  to 
complicated  configurations  and  irregular  geometrical  boundaries.  There- 
fore the  troublesome  boundary  conditions  along  the  intersection  curve  of 
intersecting  shells  is  effectively  overcome.  In  addition,  the  finite 
element  method  is  a very  convenient  and  efficient  method  for  programming 
for  electronic  computers  compared  with  other  numerical  methods  [4]. 

One  of  the  most  advantageous  applications  of  the  finite  element 
method  is  to  nonlinear  problems.  Nonlinear  behavior  can  occur  in  two 
different  forms.  The  first  is  material  nonlinearity  which  arises  because 
of  the  material  possessing  nonlinear  constitutive  laws.  The  second  is 
geometric  nonlinearity.  This  nonlinearity  is  associated  with  large 
displacements  that  cause  sufficiently  large  changes  in  the  geometry  of 
the  structure  that  the  deformed  configuration  is  used  when  writing  the 
equilibrium  conditions.  Superposition  techniques  are  no  longer  valid  for 
loadings  increased  beyond  the  proportional  limit.  However,  with  the  aid 
of  incremental  or  iterative  techniques,  the  finite  element  method  can 
handle  both  of  these  two  different  categories  of  nonlinearities  without 
major  changes  in  numerical  procedures. 

Accuracy  and  efficiency  are  two  considerations , even  essential 
issues,  that  enter  into  the  development  of  computer  programs  that  are  to 
be  applied  to  large  nonlinear  problems.  For  intersecting  cylinders,  since 
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a large  amount  of  core  storage  (or  input-output  operations  if  secondary 
devices  are  used)  and  computational  efforts  are  required  for  the  nonlinear 
solutions,  special  techniques  such  as  the  reduced  integration  concept 
should  be  considered  to  make  the  problem  tractable  in  a practical  sense. 
Very  little  mathematical  development  of  the  reduced  integration  technique 
has  been  published  to  date.  Most  publications  have  centered  around  a 
demonstration  rather  than  a development. 

1 .2  Objective  and  Scope 

It  is  the  object  of  this  study  to  develop  a general  procedure 
for  nonlinear  analysis  of  intersecting  cylinders.  The  finite  element 
method  is  selected  for  its  high  efficiency  and  convenience  in  computer 
work. 

The  progression  of  yield  is  of  particular  interest  in  this 
study.  The  three-dimensional  isoparametric  elements  are  layered  through 
the  thickness  of  the  intersecting  cylinders  in  the  region  where  the  high 
stress  gradieint  exists  while  two-dimensional  curved  shell  elements  are 
used  throughout  the  remainder  of  the  structure.  Transitional  elements 
are  employed  to  connect  these  dissimilar  three-dimensional  and  shell 
elements  together. 

Since  small  deformations  are  assumed,  only  the  material 
nonlinearity  is  considered.  The  study  is  limited  to  isotropic, 
homogeneous  materials  with  elastic,  linear  strain-hardening  behavior. 
Isotropic  strain  hardening  is  applied  with  monotonical ly  increasing 
loadings.  In  the  plastic  range,  a mixed  incremental-iterative  method 
is  included  in  the  stress  analysis.  The  Von  Mises  yield  criterion  is 
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used  to  predict  the  nonlinearity  of  elements.  The  reduced  integration 
technique  is  also  employed  to  economize  on  computer  operations. 

The  reliability  and  the  effectiveness  of  the  procedure  are 
verified  by  solving  several  examples.  Finally,  a problem  of  two  normally 
intersecting  cylinders  subjected  to  increasing  internal  pressure  is  solved. 
The  stresses  at  the  outer  and  inner  fibers  of  the  shells  are  evaluated  and 
compared  with  available  experimental  data. 

1 . j Notations 


[A] 

[B],  [B1] 

[c],  [cs] 

[D],  [D1] 

D 

d 

E 

fa 

fb 

fc 

fi 

[GU],  [GL] 
G , G 1 


H 


T3D  transformation  matrix  (Appendix  A) 

matrix  relating  nodal  displacement  and  strains,  based  on 
global  and  local  coordinate  systems 

transformation  matrix  for  3D  and  shell  transition  elements, 
respectively 

material  property  matrix  in  global  and  local  coordinates 

incremental  stress-strain  relations 

diameter  of  the  main  cylindrical  shell 

diameter  of  the  branch  pipe 

modulus  of  elasticity 

average  stress  concentration  factor 

correction  factor  when  bending  stresses  are  included 

stress  concentration  factor  of  normally  intersecting  cylinders 

body  force  components 

upper  and  lower  triangular  matrices  of  the  structural 
stiffness  matrix 

effective  area  of  the  main  cylinder  and  the  branch  pipe, 
respectively 

do/dc  ^ 
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[I] 
[0],  |J| 
J2,  J3 
[Kc] 
[K1] 
[K] 
[K$] 

[kt] 

K 

V Ma 

[N'J,  [N*] 

P 

P 

(AP}n 

{AP} 

{P} 

Po 

Q 

R 


S 


s 

ST 


unit  diagonal  matrix 

Jacobin  and  determinant  of  Jacobin 

the  second  and  the  third  invariant  of  stress  deviator 

initial  elastic  stiffness  matrix 

the  stiffness  caused  by  nonlinearity  of  material 

element  stiffness  matrix 

stiffness  matrix  of  Ahmad's  shell  element 

stiffness  matrix  of  shell  transitional  element 

element  stiffness  matrix,  and  the  hardening  parameter 

applied  moments  along  the  x',  y'  axes 

shape  functions  in  curvilinear  coordinates  for  3D  element 
and  shell  element,  respectively 

generalized  load  vector 

internal  pressure 

residual  nodal  forces  at  nth  iteration 
applied  load  increment 
load  vector 
4bo0 

distributed  surface  load 

radius  of  the  main  cylinder  also  residual  nodal  forces 
radius  of  the  branch  pipe 
stress  deviator  tensor 

nominal  hoop  stress  of  the  main  cylinder,  r/R,  and 
Ahmad's  shell  element 

nominal  hoop  stress  of  the  branch  pipe 

shell  transitional  element 

surface  traction 
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T3D 
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utop’  ubot’  umid 


V V UC 

u. 


u,  v,  w 
iT  , v' , w1 
AU 


★ 


W0 


Wo 


X,  y,  z 
x',  y',  z' 

e 

{eh  {e1 } 

{>  ♦ n » 4 

B,  a 

6W 


thickness  of  the  main  cylinder 
thickness  of  the  branch  pipe 
three-dimensional  transitional  element 
nodal  displacement  vector 

displacement  components  of  the  top,  the  bottom, 
and  the  mid-nodes  on  the  interface  of  T3D  element 

displacement  variation  along  the  edge  AB 

nodal  displacement  vector  excluding  those 
on  the  interface  of  T3D  element 

nodal  displacement  at  nodes  A,  B,  C 

departure  displacement  of  node  C (Fig.  4) 

components  of  displacement  in  the  X,  Y,  Z directions 

components  of  displacement  in  the  x',  y‘,  z'  directions 

incremental  displacement  vector 

volume  of  a given  solid  domain 

unit  direction  vector  in  the  x',  y',  z’  directions, 
respectively 

the  center  or  the  tip  deflection  when  yielding  starts 

the  center  (or  the  tip)  deflection  of  simply  supported 
beam  (or  cantilevered  beam) 

global  coordinate  system 

local  coordinate  system 

strain  vector 

strain  vector  in  global  and  local  coordinates, 
respectively 

curvilinear  coordinates 

rotations  about  the  x1,  y'  axes 

virtual  work 

stress  tensor 
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Oo  = 


a,  Aa  = 


{dee}  = 


{de^}  = 
dIP  = 


dA  = 
Y = 


P = 

y = 
[9]  = 
6 = 


a = 


3D  = 


stress  increment 

initial  yield  stress  in  simple  tensile  test 

effective  stress  and  effective  stress  increment, 
respectively 

normal  stress  of  cylindrical  shell 
elastic  strain  increment 
plastic  strain  increment 
effective  plastic  strain  increment 
nonnegative  constant 
the  rate  of  convergence 
step  length 

nondimensional  load  parameter 
Poisson's  ratio 

direction  cosine  matrix  [V-| , V ^ > V;j] 

distance  as  defined  in  Fig.  1 

linear  interpolation  factor,  intersection  angle 
of  intersecting  cylinders 

three-dimensional  element 
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CHAPTER  2 

THE  STRESS  ANALYSIS  OF  INTERSECTING  CYLINDERS 
2.1  Previous  Work 

Intersecting  cylinders  can  occur  in  a variety  of  engineering 
applications.  Therefore,  a number  of  engineering  solutions  have  been 
sought  for  these  problems  using  different  approaches  such  as  experimental, 
analytical,  or  numerical  methods.  These  procedures  have  been  performed 
to  investigate  both  the  stress  distribution  and  the  structural  behavior 
of  such  intersecting  shells.  The  previous  work  directed  to  this  problem 
is  grouped  and  summarized  below. 

2.1.1  Experimental  Work 

Experiments  conducted  on  the  intersection  region  can  commonly 
be  classed  within  two  broad  categories. 

A.  Metal  Specimens 

This  application  consists  of  measuring  the  surface  strains  at 
some  particular  points  on  an  actual  shell  or  a scale  model  machined  or 
milled  out  of  metal.  Mechanical  or  electrical  resistance  strain  gages 
are  used  for  this  purpose.  Experimental  studies  conducted  by  Mehringer 
[5]  and  Cranch  [6]  have  been  published.  The  work  presented  by  Corum  [7] 
represents  some  recent  careful  experiments.  This  latter  study  is  well 
documented  and  has  already  been  referred  to  by  other  researchers. 
Electrical  resistance  strain  gages  were  used  on  both  the  inner  and  outer 
surfaces  of  the  models  in  the  test  series.  The  series  involved  four 
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models  that  had  different  geometric  variables  and  were  subjected  to 
various  loading  patterns,  including  internal  pressure  and  end  forces 
on  the  nozzle. 

Strain  gages  were  placed  in  two  opposite  quadrants.  Each 
quadrant  had  four  lines  of  gages  which  ran  along  the  nozzle  then  when 
on  the  cylinder  radiated  from  the  junction.  The  space  between  two  gages 
on  each  line  was  based  on  the  anticipated  stress  concentration  and  on  the 
distance  from  function.  The  test  results  were  compared  with  theoretical 
predictions  derived  from  a finite  element  solution  obtained  by  using  flat 
triangular  elements. 

B.  Photoelasticity 

This  method  gives  an  overall  picture  of  stress  distribution. 

The  differences  of  principal  stresses  are  optically  measured  from  isotropic 
transparent  models  which  become  doubly  refractive  when  polarized  light  is 
passed  through  the  model.  The  newly  developed  freeze  techniques  are 
available  for  three-dimensional  models.  Upon  the  completion  of  the 
"stress  freezing"  operation,  slices  are  removed  from  the  model  and  then 
the  stresses  are  determined  by  standard  photoelastic  techniques.  4^. 

Schneider  [8]  tested  a series  of  intersecting  cylinders  which 
were  made  of  epoxy  resin  and  subjected  to  internal  pressure.  Stress 
concentration  factors,  or  stress  indices,  were  investigated  by  these 
photoelastic  tests.  Taylor  [9]  conducted  a three-dimensional  photoelastic 
study  of  stresses  around  reinforced  branch  pipe  intersections.  Taniguchi 
and  Kono  [10]  described  the  results  of  an  experimental  analysis  of  the 
nozzle  to  vessel  attachment  under  external  loadings  by  means  of  the 
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three-dimensional  photoelastic  method. 
2.1.2  Analytical  Analysis 


The  complicated  geometrical  shape  of  the  intersecting  line 
between  normally  intersecting  cylinders  creates  a difficulty  in  solving 
the  problem  analytically.  However,  the  problem  can  be  greatly  simplified 
if  the  diameter  ratio  between  the  branch  shell  and  the  main  shell  is  small. 
The  main  shell  can  be  treated  as  a shallow  shell,  so  Donnell’s  equation 
is  applicable.  The  end  section  of  the  branch  pipe  can  be  looked  upon  as 
flat.  Therefore  standard  solutions  for  cylindrical  shells  such  as  those 
presented  by  Flugge  [11]  can  be  used  directly. 

Reidelback  [12]  made  the  above  assumptions  and  derived  a 
simplified  differential  equation  to  examine  the  influence  of  internal 
pressure  on  the  elastic  behavior  of  the  intersection  region.  In  his  work, 
formulas  are  given  for  the  case  of  both  cylinders  of  equal  diameter  even 
though  the  procedure  is  valid  only  for  very  small  diameter  ratios. 

Later,  Eringen  and  Suhubi  [13]  used  Donnell's  equation  for  both 
shells  to  attack  the  same  problem,  and  established  a set  of  eight  boundary 
conditions  along  the  intersection  curve.  These  conditions  are  used  to 
determine  the  unknown  constants  of  the  analytical  solution.  The  diameter 
ratio  of  the  intersecting  cylinders  was  limited  to  less  than  one-third. 
Unfortunately,  no  numerical  example  was  presented  in  that  article. 

Bijlarrd,  Dohrmann  and  Wang  [14]  presented  results  for  the 
case  when  the  intersecting  cylinders  were  of  equal  diameter.  Thick 
shells  were  considered  and  shear  deformations  were  also  taken  into 
account.  Flugge's  equations  were  applied  to  both  cylinders  in  the 
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development  of  the  solution  process.  No  numerical  results  were  presented 
In  that  study  either. 

For  an  arbitrary  diameter  ratio  of  normally  intersecting  shells. 
Pan  and  Beckett  [15]  formulated  their  resulting  equations  on  the  basis  of 
a general  elastic  thin  shell  theory.  Donnell's  and  Flugge's  equations 
were  used  for  main  and  branch  cylinders,  respectively.  The  numerical 
example  for  the  diameter  ratio  1:2  was  selected  to  compare  with  experimental 
results.  As  pointed  out  by  Lekkerkerker  [2],  the  equations  border  on 
being  ill-conditioned  if  a numerical  procedure  such  as  collocation,  with 
points  at  equal  intervals,  is  selected  for  solving  the  equations  that 
enforce  continuity  between  the  two  shells  along  the  intersection  curve. 

Hansberry  and  Jones  [16]  also  developed  a collocation  method  to 
describe  the  elastic  behavior  of  two  normally  intersecting  cylindrical 
shells  with  small  diameter  ratios  that  are  less  than  0.2.  Their  numerical 
results  were  compared  with  the  experimental  tests  of  Cranch  and  Dally  [17]. 

2.1.3  Finite  Element  Method 

In  the  early  applications  of  the  finite  element  method  to 
intersecting  cylinders,  the  curved  shell  surfaces  were  simply  replaced 
by  flat  plate  bending  and  membrane  elements.  Because  of  discretization 
errors,  a large  number  of  such  flat  elements  was  needed  to  converge  to 
reasonable  answers. 

Prince  and  Rashid  [18]  used  triangular  plate  elements  to  solve 
the  case  of  very  thin  normally  intersecting  cylinders  with  the  diameter 
ratio  of  1:2.  Their  results  were  compared  with  experimental  data  for  a 
nozzle-to-cyl inder  intersection. 
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Greste  [19]  used  both  two-dimensional  plate  bending  and  plane 
type  elements  to  solve  the  tubular  K joint  problem. 

Hellen  and  Money  [20]  demonstrated  the  general  capabilities  of 
the  stress  analysis  program  BERSAFE  by  using  a double  layer  of  isoparametric 
elements  through  the  thickness  of  the  shells. 

Bakhrebah  and  Schnobrich  [21]  modeled  the  normally  intersecting 
cylinder  problem  by  using  three-dimensional  isoparametric  elements  along 
the  intersection  curve,  and  two-dimensional  curved  shell  elements  in  the 
regions  away  from  the  intersection.  Because  the  simple  isoparametric 
elements  formulated  by  the  displacement  method  are  inherently  too  stiff, 
incompatible  modes  [48]  and  reduced  integration  [47]  techniques  were 
investigated  as  a possible  means  for  making  the  element  and  therefore 
the  structure  more  flexible.  The  results  calculated  by  Bakhrebah  show 
good  agreement  with  the  experimental  results  obtained  by  Corum. 

The  techniques  of  nonlinear  analysis  have  been  applied  to 
structures  for  many  years.  The  application  of  employing  the  finite 
element  method  to  intersecting  shells,  however,  has  only  recently  begun. 

To  update  the  structural  stiffness  of  the  system  at  each  step  of  the 
nonlinear  analysis  is  a straightforward  but  costly  and  cumbersome 
procedure.  Some  literature  concerning  this  topic  has  been  published, 
with  several  different  approaches  being  applied. 

Mahmoud  Khojasteh-Bakht  and  Popov  [22]  provided  a general 
discussion  of  the  use  of  finite  elements  in  the  analysis  of  elastic- 
plastic  problems.  The  tangent  stiffness  method  was  employed  to  solve 
rotational  shells  subjected  to  axisymmetric  loading. 
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Gupta,  Mohraz  and  Schnobrich  [23]  used  three-dimensional 
isoparametric  elements  to  solve  a thick  circular  plate  with  circular 
openings.  Elastic-plastic  behavior  of  the  material  was  included  by 
incorporating  the  Von  Mises  yield  criterion  in  an  incremental  format. 

The  initial  stress  method  was  used  to  economize  the  evaluation  of  the 
unbalanced  nodal  forces  and  plastic  deformation  at  each  iteration. 

Larser  and  Popov  [24]  used  three-dimensional  isoparametric 
elements  for  the  elastic-plastic  analysis  of  thick-walled  pressure  vessels 
with  sharp  discontinuities  in  geometry.  A modified  incremental  method, 
termed  the  "one-step  iteration"  or  "out-of-balance  force"  method,  was  used 
to  work  out  some  numerical  examples. 

It  is  clear  from  the  above  reviews  of  the  previous  work  that  a 
reliable  general  analytical  method  for  the  nonlinear  analysis  of  cylinder- 
to-cylinder  intersections  is  not  available.  To  fill  the  need  for  an 
engineering  solution  of  the  intersecting  cylinders  problem,  the  finite 
element  method  with  its  nonlinear  feature  capable  of  representing  elastic- 
plastic  behavior  of  structures  is  highly  desirable.  Before  simulating  and 
then  discretizing  the  intersecting  cylinders  for  finite  element  models, 
some  basic  knowledge  of  the  general  behavior  of  the  structure  should  be 
known. 

2.2  Important  Geometric  Parameters  of  Intersecting  Cylinders 

The  branch  pipe  connection  is  characterized  by  the  intersection 
angle  of  the  two  cylinders,  reinforcements  around  the  intersection,  and 
three  geometric  ratios,  i.e.,  the  diameter  ratio  d/D,  the  main  vessel 
thickness  ratio  T/D,  and  the  membrane  hoop  stress  ratio  s/S  = dT/Dt. 


The  "T"  shape  connection  without  any  fillet  reinforcing  is  the  problem  of 
central  interest  in  this  study.  Its  general  behavior  is  discussed  below 
based  on  the  three  geometric  parameters. 

The  range  of  the  diameter  ratio  is  obviously  0 < d/D  £ 1. 

From  the  parametric  study  carried  out  by  Ellyin  and  Turkkan  [25],  it  was 
concluded  that  the  unreinforced  nozzle-vessel  attachment  provided  less 
strength  when  the  diameter  ratio  was  bounded  between  0.5  < d/D  < 0.6. 

The  same  conclusion  was  also  reached  by  Schroeder  [26].  For  a small  d/D 
ratio,  d/D  _<  0.2,  the  weakening  caused  by  the  cutout  in  the  main  shell  is 
relatively  small,  and  high  strength  is  anticipated.  This  has  been 
demonstrated  with  both  analytical  and  experimental  results. 

The  deformation  pattern  of  the  intersecting  cylinders  when 
subjected  to  a constant  internal  pressure  is  based  on  the  combinations  of 
thicknesses  and  radii  of  the  branch  pipe  and  the  main  shell.  If  the 
thickness  ratio  of  the  structure  is  relatively  small,  the  nominal  hoop 
stress  in  the  cylinders  is  high  but  the  distance  from  the  intersection  for 
which  the  disturbance  has  effectively  damped  out  is  small,  and  vice  versa 
for  the  large  thickness  ratio.  Therefore,  in  practical  design,  it  is 
essential  to  optimize  the  T/D  ratio  if  the  intent  is  to  use  the  material 
effectively.  The  behavior  of  thin  shells  and  that  of  thick  shells  is 
quite  different.  Accordingly,  the  analysis  approaches  are  not  the  same. 
For  thin  shells,  T/D  £ 1/20,  both  the  bending  stresses  and  the  stresses 
normal  to  the  surface  can  be  ignored,  and  only  the  membrane  stresses  due 
to  strains  in  the  middle  surface  of  the  shell  need  be  considered.  This 
is  true  except  in  the  regions  of  disturbance  such  as  penetrations, 
stiffness  changes,  and  supports.  For  thick  shells,  the  shear  effect 
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in  the  thickness  direction  is  not  negligible.  This  means  the  distortion 
across  the  thickness  invalidates  the  Kirchoff  hypothesis.  If  the  finite 
element  method  is  used  to  analyze  the  structure,  the  selection  of  element 
models  must  be  able  to  represent  the  real  behavior  of  the  shell.  In 
industrial  applications,  a T/D  ratio  in  the  range  of  1/10-1/50  is 
comparatively  common. 


2.3  Stress  Concentration  of  Intersecting  Cylinders 


The  branch  pipe  connection  consists  of  two  individual  components, 
i.e.,  the  branch  pipe  and  the  main  cylinder.  The  contact  points  of  these 
two  cylindrical  shells  form  an  intersection  curve  [16]  which,  for  the 
general  case  with  an  intersection  angle  a between  the  two  axes,  can  be 
expressed  as  (Fig.  1) 


X = r cos  <|) 

Y = r sin  <J>  - 6 sin  a 

7.  = r/I  - S2  cos2  4> 


(2.1) 


with  S = r/R  and  6 = 


R 


cos  a 


(1 


- AT 


S2  cos2  4>)  - r sin  <j>  tan  a 


where  r and  R are  the  radii  of  the  branch  pipe  and  the  main  shell, 
respectively.  A set  of  edge  forces  is  introduced  at  the  juncture  of  the 
two  shells  to  enforce  the  continuity  of  displacements  across  the 
intersection  line. 

The  standard  solution  of  the  Donell,  Flugge  or  other  form  of 
the  cylindrical  shell  equations  for  the  branch  pipe  when  subjected  to 
edge  loads  can  be  obtained  and  is  well  known.  From  such  solutions,  it  is 
evident  that  the  effect  of  an  edge  disturbance  on  a cylindrical  pressure 
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vessel  is  negligible  when 


X > 2.45  /KT 


(2.2) 


where  X is  the  distance  measured  from  the  forced  edge.  One  exception 
deserves  to  be  mentioned  here.  When  the  edge  of  the  cylinder  is  subjected 
to  a set  of  self-equilibrating  axial  forces  and  the  far  end  is  left  free, 
the  edge  disturbance  increases,  with  the  distance  away  from  the  disturbed 
edge  rather  than  dying  out.  This  causes  the  collapse  of  the  cylinder  into 
an  oval  shape.  This  phenomenon  was  first  described  by  Vlassov  [27]  in  his 
experimental  and  analytical  investigations.  Bakhrebah  also  experienced 
this  in  his  finite  element  analysis.  To  avoid  this  difficulty,  the 
constrained  boundary  conditions,  such  as  those  for  a diaphragm  closure, 
are  usually  adopted  instead  of  free  end  conditions. 

The  main  shell  is  weakened  by  the  opening  which  causes  the 
discontinuity  in  the  geometry  and  in  the  displacement  fields.  The  existence 
of  a stress  concentration  around  the  hole  can  be  visualized  by  comparing 
the  vessel  with  an  infinite  flat  plate  having  a circular  opening.  This 
plane  stress  problem  was  investigated  by  Timoshenko  [28],  It  has  been 
pointed  out  that  the  maximum  stress  is  three  times  larger  than  the  stress 
found  in  a solid  plate.  The  stress  distribution  in  the  cylinder  must  also 
be  influenced  by  the  curvature.  The  stress  concentration  varies  with  the 
size  and  shape  of  the  hole  and  may  be  three  to  four  times  as  large  as  the 
stresses  would  be  in  a solid  shell  (Taylor  [9]). 


On  the  basis  of  the  above  data,  it  is  logical  to  conclude  that 
there  is  a stress  concentration  in  the  vicinity  of  the  connection  of  a 
cyl inder-to-cylinder  intersection.  Furthermore,  experimental  data  [29] 
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have  shown  that  there  are  two  stress  peaks  along  the  intersection  line 
when  the  structure  is  under  internal  pressure: 

(1)  High  hoop  stress  in  the  vicinity  of  section  AB  due  to 
the  removal  of  material  for  the  hole  in  the  main 
cylinder  (Fig.  2). 

(2)  High  bending  stress  in  the  vicinity  of  point  C,  where 
the  internal  pressure  normal  to  the  vessel  can  only  be 
balanced  by  the  bending  action  and  the  component  from  the 
axial  force  in  the  branch  pipe.  For  large  diameter  ratios, 
the  bending  stress  increases  while  the  component  from  the 
axial  force  decreases. 

Test  results  show  [30]  that  the  bending  stress  at  point  C is 
seldom  as  high  as  the  hoop  stress  at  section  AB.  Thus  the  hoop  stress 
in  section  AB  governs  the  design  at  least  of  the  intersection.  To 
determine  the  stress  concentration  factor  at  section  AB,  an  approximate 
analysis  proposed  by  Lind  [33],  called  the  area  method,  is  available  for 
pressurized  normal  branch  pipe  connections  without  fillets  around  the 
junction.  The  area  method  avoids  rigorous  mathematic  derivations. 

Instead,  the  whole  concept  is  based  on  an  estimate  of  the  effective 
lengths  (Fig.  3)  of  the  branch  pipe  and  main  shell.  The  rate  of  decay 
of  stress  in  the  main  shell  is  assumed  to  be  a linear  variation.  The 
length  over  the  cylinder  from  the  maximum  stress  to  the  membrane  stress 
is  approximated  as  0.8  >/RT.  The  area  of  the  triangular  stress  distribution 
is  equivalent  to  the  maximum  stress  uniformly  distributed  over  an  effective 
length  0.4  /RT.  From  the  effective  lengths  of  the  branch  pipe  and  the  main 
cylinder,  the  corresponding  effective  areas  G1  and  G (effective  length  x 
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diameter.  Fig.  3)  can  be  measured  and  the  average  stress  concentration 
factor  is  computed  as 

fa  = (G/G* )/(D/2T)  (2.3) 

When  the  bending  stress  is  taken  into  account,  a correction  factor  f^  is 
introduced  as  indicated  in  Eq.  (2.4) 

fb  = 1 + (T/D)//i7S  (2.4) 

The  actual  stress  concentration  factor  f thus  is 

c 


From  a comparison  with  the  experimental  data,  the  author  quotes  a mean 
error  of  f as  less  than  3 percent  based  on  his  statistical  evaluation  of 
the  data.  Therefore,  it  is  reasonable  to  presume  that  the  high  :oncentrated 
stress  is  distributed  over  a distance  about  0.8  ^ 2.45  /Rf  from  the  junction. 

Out  of  this  region,  membrane  behavior  dominates.  This  approximation 
provides  a preliminary  estimate  for  modeling  the  structure  when  the  finite 
element  method  is  to  be  employed  to  solve  the  problem. 

2.4  Behavior  Beyond  the  Elastic  Limit 

2.4.1  Failure  Mechanism  of  Intersecting  Cylinders 

j 

The  failure  mechanism  of  intersecting  cylinders  is  essentially 
based  on  loading  conditions,  material  properties,  and  surrounding 
temperatures.  If  the  structure  is  subjected  to  a monotonically  increasing 
internal  pressure,  the  serviceability  may  end  as  a result  of  severe 
overstressing  in  some  regions.  In  other  words,  the  stress  reaches  the 
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strength  capabilities  of  the  material  in  these  regions.  But  on  the  other 
hand,  if  the  pressure  is  cyclic,  a shakedown  failure  caused  by  high  strain 
fatigue  may  be  the  controlling  factor  [31].  When  brittle  materials  are 
employed,  the  high  stress  concentration  around  the  intersection  region 
remains  right  up  to  the  breaking  point  since  the  material  has  little 
ductility  to  deform  and  hence  redistribute  more  uniformly  the  high  local 
stresses.  Therefore,  points  of  stress  concentration  along  the  intersection 
curve  have  a greater  importance  and  are  regions  of  central  interest  if 
brittle  fracture  is  a consideration.  For  ductile  materials,  a large 
deformation  may  be  developed  before  a final  plastic  rupture  occurs.  The 
environmental  conditions  also  affect  structural  behavior.  For  instance, 
the  toughness  of  intersecting  cylinders  made  from  brittle  materials  [32] 
can  be  improved  at  elevated  temperatures  even  if  the  structure  contains 
notches  or  flaws. 

Intersecting  cylinders  provide  low  serviceabilities  if  they 
undergo  little  deformation  prior  to  the  failure.  To  prevent  or  minimize 
the  brittle  fracture  possibilities  of  ductile  materials,  it  is  necessary 
to  avoid  high  stress  fields,  low  temperature  environments  and  flaws 
occurring  simultaneously.  In  most  of  the  engineering  applications  such 
as  boiler  or  nuclear  reactors,  both  the  temperature  and  the  internal 
pressure  are  very  high.  Therefore,  rupture  is  most  probably  accompanied 
by  some  plastic  deformations  if  a ductile  material  is  used.  Actually, 
from  a survey  of  the  failures  of  pressure  vessels  over  the  past  decade, 
Nichols  [34]  pointed  out  that  the  most  important  phenomenon  to  be 
considered  as  a potential  source  of  trouble  was  the  plastic  rupture  of 
the  welds,  or  ohe  weld-affected  areas  near  the  branch  attachments. 
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In  this  study,  it  is  of  particular  interest  to  investigate 
the  progression  of  plastification  after  yield  starts  but  before  the 
intersecting  cylinders  have  fractured.  To  determine  the  bounds  of 
this  range,  limit  analysis  techniques  can  be  employed  to  approximate 
the  initial  yield  load  as  well  as  the  collapse  load. 

2.4.2  Limit  Analysis 

In  limit  analysis,  the  stress-strain  relation  is  normally 
simplified  to  rigid-perfectly  plastic,  concentrating  thereby  the 
deformations  in  localized  regions.  A lower  bound  solution  is  obtained 
by  the  determination  of  a statically  admissible  system,  defined  as  any 
system  which  satisfies  the  equilibrium  conditions,  and  has  stresses  at 
every  point  at  or  below  yield.  An  upper  bound  solution  is  found  by  the 
consideration  of  a kinematically  admissible  system.  This  system  is 
defined  as  a compatible  pattern  of  displacement  for  which  the  rate  of 
external  work  is  equal  to  or  exceeds  the  rate  of  internal  dissipation. 

The  upper  bound  gives  the  maximum  collapse  load. 

The  limit  analysis  method  was  essentially  developed  for  the 
design  of  steel  frames.  Its  application  to  shells  was  first  published 
by  Drucker  [35]  in  a study  of  symmetrically  loaded  cylindrical  shells 
without  axial  forces.  Because  of  the  lack  of  rotational  symmetry, 
unlike  the  nozzle  to  spherical  shell  connections,  the  limit  analysis  of 
a branch- to-cyl indrical  vessel  is  much  more  difficult.  Attempts  at  limit 
analysis  of  this  configuration  have  only  been  made  recently.  Cloud  and 
Rodabaugh  [36]  gave  an  upper  bound  solution  for  internal  pressure  of  a 
normal  branch  pipe  connection  with  the  restriction  of  a small  diameter 
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ratio,  less  than  0.5.  Schroeder  and  Ramgarajan  [27]  also  obtained  an 
approximate  upper  bound  solution  for  diameter  ratio  between  0.4  to  1.0. 
Ellyin  and  Turkkan  [26]  have  given  a lower  bound  solution  for  internal 
pressure  by  using  the  limit  pressure  as  the  objective  function  and 
maximizing  the  objective  function  over  the  admissible  stress  field.  The 
solutions  obtained  were  over  a wide  range  of  parameters  and  were  compared 
with  experimental  data. 

To  solve  intersecting  shells  by  limit  analysis,  the  general 
geometric  relations  are  first  established.  No  simplifying  assumptions 
should  be  made  along  the  intersection  curve  if  the  solution  is  to  be  for 
the  general  case.  Then,  the  partial  differential  equations  of  equilibrium 
of  stress  resultants  are  derived.  A stress  field  which  satisfies  the 
equilibrium  equations,  all  the  boundary  conditions  and  the  stress  continuity 
condition  at  the  intersection  are  constructed  by  following  the  work  of 
Hodge  [37].  After  the  stress  fields  for  the  branch  and  main  vessel  have 
been  chosen,  a yield  criterion  is  imposed  and,  according  to  the  lower 
bound  or  upper  bound  theorem,  a set  of  inequality  conditions  are  obtained. 
The  extreme  of  the  solutions  of  these  inequality  conditions  gives  the 
lower  or  upper  bound  solution  for  the  intersecting  cylinders. 

From  a parametric  study,  it  is  found  that  the  bound  is  affected 
by  the  geometric  variables.  A relatively  small  increase  in  the  nozzle 
thickness  considerably  increases  the  limit  pressure  of  the  structure. 

This  leads  to  the  reinforcement  around  the  junction  of  cylinders  in 
practical  design.  However,  limit  analysis  gives  no  intermediate  results, 
only  the  initial  yield  load  and  the  collapse  load  of  the  branch  pipe 
connections.  To  fill  in  this  gap,  the  finite  element  approach  with 
elastic-plastic  analysis  is  desirable. 


23 


CHAPTER  3 

THE  FINITE  ELEMENT  APPROACH 


3.1  General 

The  finite  element  method  can  be  viewed  as  basically  a 
variational  approach.  When  considered  in  this  light,  the  procedures  can 
be  generalized,  thereby  extending  the  application  of  the  finite  element 
method  to  many  engineering  fields,  not  just  structural.  Using  variational 
principles,  the  governing  equations  of  a continuum  can  be  obtained  through 
the  derivation  of  a stationary  solution.  For  most  cases,  these  equations 
are  too  complex  to  be  amendable  to  closed  form  solutions  directly.  A 
usual  technique  is  to  use  the  Rayleigh-Ritz  method  to  construct  approximate 
solutions  by  reducing  or  restricting  the  unknowns  to  a small  or  finite 
number. 

The  energy  procedures  used  in  structural  mechanics  can  be 
classified  essentially  as  the  minimum  potential  energy  method  and  the 
minimum  complementary  energy  method.  The  former,  usually  referred  to  as 
the  stiffness  or  displacement  method,  associates  with  assumed  displacement 
parameters.  The  latter,  termed  the  flexibility  or  force  method,  deals 
with  a parametric  equilibrium  stress  field.  In  addition  to  these  two 
methods,  a mixed  procedure,  utilizing  the  Hell inger-Reissner  principle 
[38],  has  been  developed  by  taking  both  displacements  and  stresses  as 
primary  variables.  However,  of  these  methods  as  well  as  other  hybrid 
schemes  the  displacement  method  remains  the  most  generally  used  procedure 
in  structural  mechanics.  It  is  employed  in  the  present  study. 


The  finite  element  idealization  simulates  a real  structure  as 
an  assemblage  of  a finite  number  of  elements.  This  introduces  a 
discretization  error.  This  error  can  involve  both  the  geometry  and  the 
displacement  field.  The  upper-bounded  monotonic  convergence  is  not 
guaranteed  unless  several  sufficiency  conditions  are  satisfied.  The 
first  is  the  completeness  in  energy  that  requires  both  rigid  body  modes 
and  constant  strain  states  to  be  included  within  the  displacement  field. 
The  second  is  that  the  continuity  of  displacement  must  be  maintained 
across  any  element  interface.  However,  the  conditions  given  above  may 
be  relaxed  if  the  so-called  "patch  test"  proposed  by  Irons  [39]  is 
satisfied.  This  test  provides  a necessary  condition  for  convergence 


while  its  sufficiency  is  unproved.  Also  nothing  can  be  said  about  the 
direction  from  which  convergence  is  obtained. 


3.2  Displacement  Method 


The  matrix  formulation  of  structural  problems  arose  general 
attention  in  the  early  1950's  with  a series  of  papers  published  by  Argyris 
[40],  Turner  [41],  and  a number  of  other  investigators.  Much  progress  has 
been  made  since  then  by  introducing  new  types  of  elements  and  more 
sophisticated  computer  techniques.  Successful  developments  cover  various 
forms  of  structural  behavior  such  as  plasticity,  dynamics  and  large 
deflection  problems. 

There  are  two  basic  steps  in  the  development  of  the  displacement 


method.  One  is  the  element  formulation  and  the  other  is  the  structural 
calculation.  At  the  element  level,  the  displacements,  u,  over  the  element 
are  defined  in  terms  of  the  displacements  at  selected  points  called 
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nodal  points.  These  points  are  located  within  the  element  or  on  its 
boundary.  The  displacement  definition  is  accomplished  by  means  of 
interpolation  functions,  N: 


u = N u 


(3.1) 


where  ij  contains  all  the  nodal  displacements  of  the  element.  The  strains, 

£,  at  any  point  in  the  element  are  obtained  by  taking  appropriate  derivatives 
of  displacement  field  with  respect  to  the  selected  element  coordinate  system. 
The  strain-displacement  relations  can  be  expressed  as 


e = B u 


(3.2) 


where  the  coefficients  of  matrix  B are  functions  of  the  nodal  coordinates. 

The  condition  of  equilibrium  is  obtained  by  applying  the  principle 
of  minimum  potential  energy 


6(W.  + W. ) =0 
l e 


(3.3) 


where  can  be  expressed  as  the  integral  of  strain  energy  over  the  voli 


of  element  under  deformation 


Wi  = | 2 Di jk£  ei j 


ek*  dV 


\ uT  BT  DB  u dV 


(3.4) 


The  external  work  done  by  the  surface  traction  T.  and  the  body  force  f^ 


is  given  as 


We  - u1  dV  - u.  ds 


(3.5) 


Substituting  these  into  Eq.  (3.3)  and  taking  the  first  variation,  the 
following  equation  is  obtained 


P = NT  f dV  + NT  T ds  (3.8) 


Jv  Js 

Since  the  virtual  displacements  6u^  are  arbitrary,  Eq  (3.6)  can  be 
simplified  as 

K u = P (3.9) 


K and  P are  the  stiffness  matrix  and  the  generalized  load  vector  of  the 
element,  respectively. 

At  the  structural  level,  the  total  structural  stiffness  matrix 
and  the  structural  load  vector  are  set  up  following  the  superposition 
technique  to  assemble  all  the  elements  of  structure  together  properly. 
After  incorporating  the  boundary  conditions,  the  nodal  displacements  can 
be  solved.  The  strains  and  the  stresses  can  therefore  be  evaluated 
wherever  desired. 

The  structural  stiffness  matrix  is  characterized  by  being 
symmetric,  banded,  sparsely  populated  and  positive  semi-definite.  Only 
the  upper  or  the  lower  triangular  form  obtained  by  decomposition  is 
considered  in  computation  with  only  a half  band  of  it  being  stored. 

3.3  Structural  Modeling 


stress  gradient,  a fine  mesh  is  needed  to  achieve  accuracy  in  results. 
A coarse  mesh  may  cause  the  violation  of  local  equilibrium  even  at 
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integration  points.  These  are  the  points  that  are  used  as  control  points 
to  predict  the  occurrence  of  the  plastic  actions  in  the  structure  when 
nonlinear  analysis  is  being  used.  This  drift  from  the  true  response  can 
be  kept  to  a minimum  if  the  structure  is  modeled  properly. 

Hand,  et  al.  [42]  used  a layered  concept  to  investigate  the 
progression  of  cracking  that  develops  in  a concrete  slab  or  shell  under 
external  loads.  The  thickness  of  each  element  is  divided  into  several 
layers.  Each  layer,  in  turn,  may  have  different  material  properties. 

The  nodal  displacements  are  converted  to  middle  surface  strains  and 
curvatures,  then  to  layer  strains  by  employing  the  Kirchoff  assumption 
that  implies  normals  to  the  middle  surface  remain  straight  and  normal 
after  deformation.  From  the  stress  calculation,  the  excess  stresses  in 
each  layer  are  accumulated  and  converted  back  as  unbalanced  nodal  forces 
for  the  next  iteration.  This  procedure  worked  well  for  the  plate  and 
smooth  shallow  shells  that  were  studied.  A layered  concept  which  allows 
the  plasticity  to  propagate  through  the  thickness  as  well  as  along  the 
surface  will  be  employed. 

The  intersecting  shell  is  the  problem  of  particular  interest  in 
the  present  study.  From  previous  studies,  it  is  known  that  the  stresses 
in  intersecting  shells  decay  sharply  to  reach  the  membrane  stress  levels 
away  from  the  intersection  region.  Therefore,  a layered  system  need  to 
be  considered  only  in  the  vicinity  of  the  intersection.  The  desirability 
of  restricting  the  layering  to  as  small  an  area  as  possible  is  to  economize 
the  computational  effort.  Because  of  the  complicated  geometrical  shape  of 
the  intersection  curve  and  the  displacement  variation  through  the  thickness 
in  this  region,  it  is  undesirable  to  impose  the  Kirchoff's  assumption  on 
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displacement  field  as  Hand  did.  The  three-dimensional  isoparametric 
element  family  is  therefore  used  to  divide  the  thickness  of  the  structure 
into  several  layers  in  this  intersection  region. 

To  avoid  an  abrupt  change  of  stiffness  in  the  structure  away 
from  this  region,  a method  of  grading  the  mesh  from  a fine  to  a coarse 
mesh  is  obtained  by  connecting  every  two  layers  of  the  previous  grid  to 
a single  layer  in  the  adjacent  grid.  The  same  procedures  are  repeated 
until  only  one  layer  represents  the  entire  thickness.  Then  the  shell 
transitional  elements  are  used  to  connect  to  the  two-dimensional  curved 
shell  elements  which  are  used  throughout  the  remainder  of  the  structure. 
Although  the  modeling  method  is  developed  for  interesecting  shell  analysis, 
the  general  nature  of  the  procedures  used  is  applicable  to  other  kinds  of 
structures  as  wel 1 . 

3. 4 Element  Stiffness  Matrix 

3.4.1  The  Isoparametric  Family 

Three-dimensional  solid  elements  are  capable  of  correctly 
representing  the  behavior  of  a beam,  plate,  or  shell  including  any  of 
the  varied  aspects  of  structural  components  because  they  enable  bodies 
with  curved  boundaries  to  be  treated  with  a limited  number  of  elements. 

A general  isoparametric  element  suggested  by  Irons  [43]  is  adopted  in 
the  present  study.  With  that  formulation  it  is  possible  to  add  any 
number  of  intermediate  nodes  to  the  individual  edges  of  an  eight-node 
brick  element  by  employing  the  so-called  "departure  concept".  The 
displacement  variables  at  intermediate  nodes  are  treated  as  the 
difference  or  departure  from  the  linear  displacement  variation 
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between  two  corner  nodes. 

One  way  of  expressing  the  displacements  at  the  edge  of  an 
element  would  be  to  use  interpolation  functions  in  the  form  of  a quadratic 
variation  on  the  edge  AB  (Fig.  4).  This  can  be  expressed  as 

uAB  • Va  + Vb  + Ncuc  (3J0) 

where 

AB 

U = quadratic  response  along  the  edge  AB 
Ua>  Ug,  Uc  = nodal  values  at  nodes  A,  B and  C,  respectively 
Na,  Ng  = quadratic  shape  functions  corresponding  to  nodes 
A,  B and  C,  respectively 

On  the  other  hand,  the  displacements  of  the  intermediate  nodes  can  be 
written  as 

UC  " °C  * J<.UA  ♦ V (3-1” 

where  is  the  departure  displacement  of  node  C as  shown  in  Fig.  4. 

Substitution  of  Eq.  (3.11)  into  Eq.  (3.10)  gives  an  alternate 
way  of  expressing  the  displacements  as 

uAB  * Wa  * ¥b  + ¥c  (3J2) 

where 

*A  = NA  + \ NC  = I « 

^B  = NB  + 2 NC  = \ ^ + ^ 

nc  = (i  - e) 
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and  Ng  are  the  linear  shape  functions  of  the  corner  nodes  A and  B 
(C  = ±1).  With  or  without  the  last  term  of  Eq.  (3.12),  the  variations 
along  the  edge  AB  become  quadratic  or  linear.  This  means  whenever  an 
additional  intermediate  node  is  introduced  to  any  edge  of  an  8-node 
element,  only  the  last  term  needs  to  be  added  without  changing  the  rest 
of  the  equation.  In  a manner  similar  to  Eq.  (3.12),  more  intermediate 
nodes  can  be  incorporated  to  the  edge  AB  in  order  to  define  a higher 
degree  of  response. 

The  isoparametric  displacement  field  within  an  element  is 

given  as 


r ^ 
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(3.13) 


i i 


\ t 

n 


t* 


•* 


I 


For  a general  curvilinear  element  the  geometric  transformation  relationship 
between  the  global  cartesian  coordinates  and  the  local  isoparametric 
coordinates  is  established  by  Eq.  (3.14)  as 


(3.14) 


y)=  l M*i 


The  strain-displacement  relation  is  defined  by  proper  differentiation  of 
the  displacement  field  as 
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(3.15) 
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Substituting  into  Eq. 


[K] 


(3.7),  the  element  stiffness  is  obtained 
[B]T[D][B]  dV  (3.16) 


where  [D]  is  the  material  property  matrix  defined  as  the  stress-strain 
relationship  of  a homogeneous  linearly  elastic  material.  This  matrix  is 
given  in  Appendix  A.  The  volume  element  dV  has  to  be  transformed  to  a 
curvilinear  coordinate  system  for  the  integration  process 


dV  = dx  dy  dz  = |J|  d£  dn  dc 


(3.17) 


where  |J|  is  the  determinant  of  the  Jacobian  matrix.  Equation  (3.16)  is 
now  of  the  form 


[K] 


fl 

rl 

-1  • 

-1  • 

[B]T[D][B]  | J | d£  dn  dc 


(3.18) 


Both  [B]  and  [D]  contain  many  null  factors.  A lot  of  intermediate 
calculations  can  be  eliminated  if  the  calculation  of  the  element  stiffness 
matrix  is  broken  into  parts  and  only  the  non-zero  terms  are  executed.  A 
more  detailed  discussion  of  this  is  present  in  Appendix  A. 

There  are  several  points  that  bear  mentioning  here. 

1.  The  requirements  of  continuity  and  those  for  the  constant 
strain  states  are  satisfied  in  the  isoparametric  element 
family,  thus  insuring  convergence. 

2.  The  elements,  however,  are  far  too  stiff  against  flexure, 
especially  the  low  order  elements.  For  instance,  a three- 
dimensional  8-node  brick  type  isoparametric  element  (Fig.  5) 
develops  the  constraints  on  the  transverse  displacement 
mode  because  of  the  appearance  of  parasitic  shear  strain 
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energy.  Therefore,  softening  procedures  such  as  the  use 
of  nonconforming  modes  or  a reduced  integration  technique 
are  usually  utilized  to  overcome  this  deficiency. 

3.  The  Kirchoff's  hypothesis,  i.e.,  the  normal  to  the  mid- 
surface remains  straight  and  normal  after  deformation, 
is  adequate  for  most  shell  problems.  Thus,  the  unnecessary 
nodes  through  the  thickness  of  an  element  can  be  eliminated 
to  minimize  computational  efforts.  In  addition,  since  the 
concept  of  a layered  system  is  employed,  the  high  order 
displacement  variation  along  the  shell  thickness  direction 
can  be  approximated  by  several  elements  having  only  two 
nodes  through  the  element  thickness. 

3.4.2  Three-Dimensional  Transitional  Element 


It  is  desired  to  connect  both  elements  A and  8 to  only  a single 
element  in  the  next  region  (Fig.  6).  Obviously,  the  continuity  across  the 
element  interface  will  be  violated  if  a quadratic  displacement  variation 
is  allowed  in  element  C while  in  elements  A and  B the  displacement 
variation  on  the  corresponding  face  is  only  linear.  Therefore,  it  is 
reasonable  to  impose  a linear  displacement  variation  on  the  interface 
between  these  three  elements.  With  this,  the  constraint  equation  of  the 
mid-nodes  can  be  established  as 


{u} 


mid 


({u} 


top 


lulbot> 


(3.19) 


where  [C]  is  a transformation  matrix  and  is  described  in  Appendix  A. 
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Introducing  this  relationship  into  the  formulation  of  the  [B] 
matrix,  we  get 

(3.20) 


(3.21) 

In  this,  [I]  is  the  unit  diagonal  submatrix  and  U is  the  nodal  displacement 
vector  excluding  those  on  the  interface  of  T3D  element. 

The  element  stiffness  matrix,  independent  of  nodes  4,  5,  and  6 
(Fig.  6),  for  elements  A and  B can  be  obtained  through  the  substitution  of 
[B‘]  for  [B]  in  Eq.  (3.18).  With  this  transition  element  it  is  now  possible 
to  connect  the  layered  region  to  any  general  three-dimensional  element 
with  3 nodes  on  both  the  top  and  the  bottom  edges  of  the  interface. 

3.4.3  Ahmad's  Shell  Element 

Ahmad's  shell  element  [44]  is  extracted  from  a 16-node  three- 
dimensional  isoparametric  element.  Conversion  to  the  shell  element 
precludes  the  possible  ill  conditioning  that  occurs  when  the  shell 
thickness  is  very  small  compared  to  the  other  dimensions  of  the  element. 
Furthermore,  it  reduces  the  number  of  unknowns.  In  this  element,  the 
constraint  of  straight  normals  is  imposed  and  the  strain  energy 
corresponding  to  normal  stresses  perpendicular  to  the  middle  surface 
is  ignored. 
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For  a typical  thick  shell  element  (Fig.  7),  the  geometric  shape 
can  be  defined  as  Eq.  (3.22)  if  there  are  no  intermediate  nodes  in  the 


element  thickness  direction. 
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where  N^U,  n)  are  the  shape  functions  for  a two-dimensional  quadrilateral 

element  as  shown  in  Fig.  8.  For  convenience,  Eq.  (3.22)  can  also  be 

rewritten  in  a form  specified  by  a nodal  vector  that  connects  the  pairs  of 

nodes  i.  „ and  i.  ..  and  the  midsurface  coordinates 
top  bottom 
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Similarly,  the  displacement  field  at  any  point  in  the  element 


can  be  expressed  as 
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where  Aui , Av^  and  Awi  are  relative  displacements  or  displacement 
differences  of  i'top  and  ibottoni-  Awi  is  ignored  as  is  the  corresponding 
strain  energy  in  the  thickness  direction.  Equation  (3.25)  can  therefore 
be  rewritten  as  (Fig.  9) 
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(3.26) 


mid 


or 


v > = 


(3.27) 


wi  th 


<v  = 


VT 

wi  } 
a. 

vB, 


i > 


mid 


[N*]^,  expressed  explicitly  in  Appendix  B.  u^. , vi  and  w^  are 
midsurface  nodal  displacements  along  the  three  global  coordinate  axes 
while  ai  and  B.  are  rotations  about  V2i  and  ^ . , respectively.  In  this. 


V1i  ' 1 x V3i 


V2f  * hi  * hi 


(3.28) 


where  i is  a unit  vector  along  the  global  x axis.  By  this,  three  local 
Cartesian  coordinate  axes  x',  y'  and  z'  are  defined  at  the  midsurface 
node  i . 
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The  strain  components  of  interest  at  any  point  in  the  element 


are  established  based  on  its  local  coordinate  system. 
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(3.29) 


After  the  [D1]  matrix  of  an  anisotropic  material  has  been 
constructed  (Appendix  B)  the  stiffness  matrix  can  be  found  in  a systematic 
manner  following  some  coordinate  transformations. 
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-1  J-l  J-l 


C B ' ]T[D' ][B ' ] | J | d?  dn  dc 


(3.30) 


A more  detailed  description  of  the  formulation  of  the  [B‘]  matrix  is 
given  in  Appendix  B. 

With  the  Ahmad  element,  there  is  the  inherent  weakness  of  being 
far  too  stiff  against  bending.  This  weakness  stems  from  an  excessive 
extraneous  shear  strain  energy  which,  in  turn,  causes  much  slower 
convergence  than  desired.  As  pointed  out  by  Pawsey  [45],  this  weakness 
can  be  diminished  by  the  use  of  the  reduced  integration  method.  When 
reduced  integration  is  used,  the  stiffness  matrix  may  become  singular  as 
zero  strain  energy  may  occur  at  integration  points.  This  is  for  the 
individual  elements.  But  when  assembled  into  a structure,  the  singularity 
is  suppressed  upon  joining  the  element  to  others. 


37 


3.4.4  Shell  Transitional  Element 

In  selecting  the  structural  modeling  to  be  used,  it  becomes 
desirable  to  connect  together  two  dissimilar  elements,  the  three- 
dimensional  solid  element  and  the  two-dimensional  curved  shell  element. 
The  transition  element  provides  this  service.  It  converts  the  five 
shell  degrees-of-freedom  at  the  edge  where  it  is  intended  to  connect  to 
a three-dimensional  element  back  to  six  degrees-of-freedom,  three  each  at 
the  top  and  the  bottom  face  of  the  element. 

Since  the  departure  concept  is  employed  with  the  three- 
dimensional  element,  the  same  modification  has  to  be  applied  to  the 
transition  element  at  each  node  that  is  intended  to  connect  to  a three- 


dimensional element  node  for  which  the  departure  is  used.  The  other 
nodes  remain  the  same  as  in  the  Ahmad  shell  element.  After  the  new  shape 
functions  for  these  nodes  have  been  defined,  as  described  in  Appendix  B, 
a new  shell  element  stiffness  matrix  can  be  generated  by  standard 
procedures  as  in  the  Ahmad  element.  Then,  the  stiffness  matrix  of  a 
transition  element  is  obtained  from  Eq.  (3.31). 

[Kt]  = [C]T[Ks][C]  (3.31) 

where  [C]  is  defined  as 

{Us>  = [C]{u3D)  (3.32) 

The  nodal  forces  that  correspond  to  the  displacement  { Uj}  in 
transition  elements  are  obtained  by 
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The  derivation  of  matrix  [C]  is  also  described  in  Appendix  B. 

3.5  Generalized  Loads 

In  the  finite  element  system,  loads  are  prescribed  only  at  the 
nodal  points  and  in  the  directions  corresponding  to  displacement  components. 
When  distributed  loads  are  applied  to  the  structure,  the  equivalent 
generalized  nodal  loads  as  outlined  in  Section  3.2  are  used  because  of 
their  computational  efficiency.  The  distributed  pressure  load  is  of 
particular  interest  and  is  described  here. 

3.5.1  Three-Dimensional  Element 


The  pressure  load  Q is  applied  to  either  the  top  or  the  bottom 
face.  Let  ri  be  the  unit  vector  normal  to  the  surface  at  any  point  p. 

The  equivalent  nodal  loads  can  be  established  on  the  basis  of  the 
equivalence  of  the  work  done  during  a virtual  displacement  consistent 
with  that  for  the  distributed  load  Q. 


or 


<5W  = {6ui}T(Pi> 


6 { U }T  Q n dS 
s 


= {<$u.j} 


[N]T  Q n dS 
s 


{V 


n dS 


(3.34) 


3.5.2  Ahmad's  Element 


The  same  procedures  are  followed  for  the  shell  element.  The 
generalized  load  vector,  at  node  i,  of  a pressure  load  Q on  the  surface 
z,  = 1 in  Ahmad  element  can  be  constructed  as 


3.6  Reduced  Integration  Technique 

Obtaining  the  stiffness  matrix  of  an  element  involves  an 
integration  over  the  volume  of  the  element  as  Eq.  (3.18)  indicates. 

For  most  cases,  the  form  of  these  integrals  is  far  too  complex  to  be 
carried  out  explicitly.  To  circumvent  this,  these  integrations  are 
frequently  done  numerically.  The  quadrature  rule,  as  the  sume  of  a 
series  of  products  of  weighting  coefficients  times  the  value  of  the 
integrand  evaluated  at  a number  of  points  is  used.  Obviously,  the  fewer 
the  number  of  points  involved,  the  less  the  amount  of  computation  required. 
From  numerical  experimentation,  it  has  been  shown  that  less  accurate 
numerical  integration  rules  can  produce  better  displacement  and  stress 
values.  This  happens  because  of  introducing  compensating  errors.  The 
reduced  integration  technique  was  first  employed  by  Doherty,  et  al.  [46] 
on  plane  quadrilateral  elements  and  later  by  Pawsey  [44]  on  curved 
elements.  Zienkiewicz,  Taylor  and  Tuo  [47]  used  a general  reduction  in 
Gaussian  integration  order  rather  than  applying  the  reduction  only  on 
the  shearing  strain  energy  components.  They  demonstrate  much  better 
accuracy  than  in  Ahmad's  first  work.  Choi  and  Schnobrich  [48]  compared 
the  results  obtained  by  including  nonconforming  modes  with  those  by  the 
reduced  integration  technique.  Dovey  [49]  investigated  the  applicability 


of  three-dimensional  elements  for  general  shell  problems  by  employing 
reduced  integration  procedures. 

Theoretically,  when  the  limiting  subdivision  of  a structure 
is  approached,  each  element  approaches  a state  of  constant  energy.  The 
total  energy  is  then  obtained  by  summing  these  constant  energies  over 
all  the  element  volumes.  Therefore,  to  obtain  a minimal  degree  of 
accuracy,  the  quadrature  rule  must  be  able  to  evaluate  the  element  volume 
correctly.  For  three-dimensional  isoparametric  elements,  a second  order 
quadrature  rule  meets  this  requirement  and  converges  rapidly.  The 
improvement  of  element  performance,  when  using  a 2 by  2 integration  rule, 
is  attributed  primarily  to  the  elimination  of  the  extraneous  shear  strain 
energy  at  the  ordinates  of  the  two  Gaussian  integration  points.  That 
makes  the  element  far  from  being  too  stiff. 

In  the  finite  element  idealization,  a geometric  regularity 
condition  is  imposed  on  the  element  in  a practical  mesh  subdivision 
process.  Even  for  quite  irregular  configurations,  numerical  examples 
show  the  volume  error  induced  is  very  small  compared  with  the  other 
approximations  involved.  Thus,  any  order  of  quadrature  rule  actually 
yields  the  correct  volume  as  the  element  limit  is  reached.  This  insures 
the  convergence.  A reversed  argument  was  therefore  proposed  by  Dovey 
that  the  convergence  might  be  assured  if  a positive  semidefinite  stiffness 
matrix  of  appropriate  rank  was  obtained.  In  other  words,  any  reduced 
integration  scheme  considered  should  be  able  to  insure  convergence  as 
long  as  it  results  in  a stiffness  matrix  that  is  positive  definite  after 
the  rigid  body  modes  have  been  removed.  No  rigorous  analysis  was 
provided  but  rather  justification  was  based  on  numerical  examples. 
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However,  a reduced  integration  method  tends  to  reduce  the  value  of 
stiffness  of  an  element  to  below  the  value  evaluated  exactly.  This 
softening  allows  the  use  of  a coarser  mesh  and  economizes  the  cost.  The 
reduced  integration  technique  is  now  widely  used  in  research  work  even 
though  the  merit  of  the  monotonic  convergence  property  is  lost.  In  the 
present  study,  a second  order  integration  is  utilized  for  all  the  elements. 

3.7  Equation  Solver 

After  the  structural  stiffness  matrix  has  been  generated,  it  is 
desirable  to  solve  for  the  nodal  displacements.  There  are  several 
numerical  schemes  available  for  this  purpose.  The  Gaussian  elimination 
method  is  often  used  for  large  systems  due  to  its  efficiency,  and  it  is 
adopted  in  the  present  study.  In  the  method,  the  whole  process  is  divided 
info  three  separate  steps,  i.e.,  the  decomposition,  the  forward  substitution, 
and  the  backward  substitution. 

In  the  decomposition,  the  structural  stiffness  matrix  is  split 
into  upper  and  lower  triangular  matrices  as 

[K]  = [GL][GU]  (3.36) 

where  [G^1]  is  nearly  the  transpose  of  [ G*- ] except  for  the  fact  that  [G^1] 
has  been  normalized  to  make  all  diagonal  terms  equal  to  unity.  Because 
of  a large  amount  of  zero  factors  scattered  inside  the  banded  structural 
stiffness  matrix,  the  efficiency  can  be  improved  by  bookkeeping  the  first 
nonzero  entry  in  each  row  to  avoid  unnecessary  computer  operations. 

Because  of  the  symmetry  of  the  structural  stiffness  matrix,  only  the 
lower  triangle  is  developed  and  stored  in  blocks  on  the  secondary  devices 
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of  the  computer  system.  Each  block,  with  the  same  size  of  nodal  degree- 
of-freedom,  is  brought  back  to  the  main  core  in  turn  and  broken  into 
submatrices  for  execution.  Inversions  of  the  diagonal  submatrices  are 
necessary  and  the  singularity  caused  by  any  mistake  is  detected. 

The  forward  substitution  computes  the  intermediate  results  {x}. 

[GL]{x>  = {P}  (3.37) 

The  backward  substitution  gives  the  nodal  displacement  by  operating 

[ GU ] { U } = ( x}  (3.38) 

Only  the  forward  and  the  backward  substitutions  need  to  be  carried  out 
when  the  solutions  of  different  load  vectors  are  required.  This  feature 
provides  the  feasibility  for  nonlinear  analysis  when  the  iteration  method 
is  employed. 


CHAPTER  4 


PLASTIC  ANALYSIS 


4.1  General 

Many  engineering  applications  require  yielding  as  an  essential 
phenomenon  of  a structural  material.  Beyond  the  yield  point,  the  load- 
displacement  relation  is  no  longer  proportional. 

Plastic  deformation  is  the  movement  of  one  layer  of  atoms  with 
respect  to  another  inside  the  material.  This  slipping  is  often  associated 
with  the  presence  of  shearing  forces.  From  a series  of  tests  with  ductile 
materials  conducted  by  Bridgman  [50],  it  has  been  concluded  that  yielding 
does  not  occur  under  hydrostatic  stress  even  though  those  stresses  may  be 
very  high.  For  the  hydrostatic  state,  no  shear  stress  occurs  in  any 
direction.  With  this  experimental  observation,  the  mathematical  models 
for  plasticity  can  be  considerably  simplified. 

Plastic  deformations  are  irrecoverable.  Also,  the  strain, 
unlike  that  considered  in  elasticity,  is  not  uniquely  determined  by  the 
final  stress  but  depends  instead  on  the  loading  path.  Incremental  theory 
[51]  is  thus  necessary  to  relate  strain  increments  to  stresses.  Deforma- 
tion theory,  which  determines  the  total  strain  components  in  terms  of  the 
state  of  stress,  will  not  be  considered  here. 

Lack  of  strain  recovery  is  caused  by  "locked  in"  residual 
microstresses  which  result  in  a Bauschinger  effect  upon  unloading.  For 
an  isotropic  hardening  material,  however,  the  Bauschinger  effect  is 
ignored. 
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4.2  Yield  Criteria  and  Incremental  Theory 

For  a simple  member  with  uniaxial  load,  yielding  is  easy  to 
determine.  But  for  most  structures,  yielding  may  be  caused  by  a 
combination  of  stress  components  a...  On  the  basis  of  experimental  work 

' J 

and  theoretical  modifications,  numerous  yield  criteria  have  been  proposed 
to  describe  the  behavior  of  material  after  yielding  has  occurred. 

In  general,  the  yield  criterion  depends  upon  the  state  of 
stress  at  the  point  under  consideration.  Therefore,  the  condition  that 
a material  has  been  loaded  to  the  initial  yield  can  be  expressed  as 

F(ai j)  = K (4.1) 

where  F = the  loading  function,  and 

K = the  hardening  parameter  which  describes  the  strain  history. 

This  equation  represents  a yield  surface  in  six-dimensional  stress  space. 
Some  materials,  such  as  metals  or  crystalline  rocks,  yield  no  plastic 
volume  change  during  plastification.  The  hydrostatic  or  spherical  stress 
state,  as  experiments  show,  does  not  cause  any  plastic  deformation. 

Hence,  it  is  usual  to  substract  the  hydrostatic  component  from  actual 
stresses  and  use  only  the  remaining  stress  deviator  in  the  yield  function. 
Equation  (4.1)  can  thus  be  written  in  a general  form  as 

F(J2,  J3)  = K (4.2) 

where  and  are  the  invariants  of  the  stress  deviators. 

The  Von  Mises  [51]  and  the  Tresca  [51]  yield  criteria  are  the 
most  widely  used  for  ductile  material.  However,  there  are  several 
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drawbacks  with  the  Tresca  criterion.  First  of  all,  the  principal  stresses 
have  to  be  known.  Otherwise,  the  Tresca  yield  function  becomes  quite 
complicated  compared  to  the  Von  Mises  criterion.  Secondly,  with  Tresca, 
the  intermediate  principal  stress  has  no  effect  on  yielding  while  in  the 
Von  Mises  criterion  all  three  principal  stresses  are  taken  into  account. 
From  experimental  data  plots,  it  has  been  found  that  the  Von  Mises 
criterion  generally  provides  closer  correlation.  Thirdly,  in  order  to 
find  the  maximum  shear  stress  at  one  point,  it  is  necessary  to  make 
comparisons  continuously  in  order  to  find  the  order  of  principal  stresses. 
This  makes  the  use  of  the  Tresca  procedure  less  desirable.  In  the  present 
study  therefore,  the  Von  Mises  criterion  is  used  to  investigate  the 
plasticity  of  the  shell  material. 

For  the  Von  Mises  criterion,  Eq.  (4.2)  can  be  expressed  in  a 
simple  form  as 


,12 

J2  3 a° 


(4.3) 


where  a0  is  the  initial  yield  stress  in  simple  tension.  A more  detailed 
derivation  of  Eq.  (4.3)  is  shown  in  Appendix  C. 

The  yield  surface  of  a Von  Mises  criterion  in  the  stress  space 
can  be  interpreted  geometrically  as  a circular  cylinder  with  its  axis 
equally  inclined  to  the  stress  axes.  For  an  isotropic  hardening  material, 
the  yield  cylinder  expands  without  changing  its  shape  as  the  loading  is 
increased.  This  can  be  expressed  as  in  Eq.  (4.4). 


dF  = -g-  do..  > 0 
8oij  1J  " 


(4.4) 


If  dF  equals  zero,  the  equation  indicates  a neutral  loading  case. 
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This  means  the  stress  state  is  moving  on  the  yield  surface.  The  condition, 
dF  > 0,  means  the  controlling  stress  state,  in  other  words  the  yield 
surface,  is  expanding. 

For  subsequent  yield  surfaces,  after  the  initial  yield  has 
occurred,  it  is  operationally  desirable  to  correlate  the  state  of  stress 
with  the  history  of  deformation  by  a single  curve  as  in  the  result  of  a 
simple  tensile  test.  The  effective  stress  and  the  effective  plastic  strain 
increments  are  thus  introduced.  Their  definition  is  as  expressed  in  Eq. 
(4.5)  and  Eq.  (4.6),  respectively. 


0 =4 (sio  si/  <4-5> 

dIP=/|  (dejj  d£^)%  (4.6) 


where  S..  is  the  stress  deviator  tensor.  Differentiating  both  sides  of 
Eq.  (4.5),  we  obtain 


do  = {f^-}  {do} 
da 


(4.7) 


To  obtain  general  stress-strain  relations  after  yielding  has 
started,  the  Prandtl -Reuss  assumption  [51]  is  employed.  This  theory 

p 

states  that  plastic  strain  increments  {de  } are  proportional  to  the 
instantaneous  stress  deviation,  i.e., 

{deP}  = S..  dX  = {|2. } d?P  (4.8) 


The  stress  increments  {de}  are  now  related  to  the  elastic  strain  increments 
{deg}  through  Hooke's  law 

{do}  = [D]{deg}  = [D]({de}  - {deP}) 


(4.9) 
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Substituting  Eq.  (4.8)  into  Eq.  (4.9)  and  premultiplying  both  sides  of 

an  T 

Eq.  (4.9)  by  {|^-}  as  Zienkiewicz  [52]  did,  we  get 

da  = {|§}T[D]({de}  - de  P)  (4 


<ff> 


{ de } = [W]{dc} 


where  H = and  can  be  obtained  for  a given  stress-strain  curve.  With 
de 

the  substitution  of  Eqs.  (4.8)  and  (4.11)  into  Eq.  (4.9),  it  can  be 


rewritten  as 


{da}  = ([D]  - [D]{f|}[W])  {de}  = [Dep]{de} 


[0  ] varies  with  the  state  of  stress  and/or  the  deformations. 

Therefore,  the  structural  stiffness  corresponding  to  the  current  material 
properties  becomes  a function  of  the  existing  displacements,  i.e., 

[ K ( U ) ] { AU } = {AP}  (4.13) 

This  means  the  load-displacement  relations  are  nonlinear.  When  loads  are 
applied,  the  calculated  load-displacement  relation  departs  from  the  proper 
curve  and  corrective  procedures  have  to  be  employed  in  order  to  bring  the 
solutions  back  to  satisfying  the  equilibrium  equations. 

4 . 3 Solution  Method 

The  problem  defined  by  Eq.  (4.13)  can  also  be  considered  as  a 
system  of  n simultaneous  independent  nonlinear  equations  of  the  form  [53], 
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where 


i = l ,n 


= 0 


c "\ 


(4.14) 


{ u ^ } are  components  in  a n-dimensional  space  or  nodal  parameters  in  the 
finite  element  method.  The  norm  of  the  vector  in  Eq.  (4.14)  is 

F(U)  = ||  MU)  ||2  = ? MU)  >0  (4.15) 

1 i=l  1 

Equation  (4.15)  becomes  zero  only  if  each  f^(U)  = 0,  and  this  can  happen 
only  if  U satisfies  Eq.  (4.14).  Hence,  the  problem  of  finding  a nonlinear 
solution  is  equivalent  to  the  problem  of  searching  for  the  minimum  of  Eq. 
(4.15).  By  employing  a truncated  Taylor's  series  expansion,  Eq.  (4.15) 
can  be  expressed  as  a quadratic  approximation  in  terms  of  AU. 

F(U)  = F(U0)  + VTF(U0)  AU  + \ f T V2F(U0)  AU  (4.16) 


Differentiating  F(U)  with  respect  to  each  of  the  components  of  AU  and 
equating  the  resulting  expression  to  zero,  we  get 

AU  = [V2F(Uo)]_1  [VTF(Uo)]  (4.17a) 

where  [V2F(U0)]_1  is  the  inverse  of  the  Hessian  matrix  defined  as  the 
matrix  of  the  second  partial  derivatives  of  F(U)  with  respect  to  U 
evaluated  at  U0.  The  inverse  of  the  matrix  should  be  positive  definite 
for  all  AU  i 0.  The  application  of  Eq.  (4.17a) in  the  finite  element 
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method  reverts  to  an  equation  of  the  form 

AU  = [K  + K']-1  {AQ}  (4.17b) 

where  [K1]  is  the  stiffness  caused  by  nonlinearity  and  {AQ}  are  the 
incremental  nodal  forces.  [Kc]  is  the  initial  elastic  stiffness. 

The  solution  process  indicated  by  Eq.  (4.17)  is  the  Newton- 
Raphson  method  (64).  However,  there  is  a serious  drawback  to  this 
method.  The  inversion  of  [K]  at  each  iteration  makes  the  procedure 
very  inefficient  for  large  systems.  In  the  modified  Newton-Raphson 
method  [54],  the  continual  requirement  for  carrying  out  the  inversion 
is  avoided  by  using  the  same  stiffness  throughout;  i.e.. 


AU  = LKc]-1  {AQ} 

(4.18) 

The  new  U is  updated 

U]  = U0  + AU 

(4.19) 

or 

U]  = U0  + A0N 

(4.20) 

where 

A0  = ||  AUo  ||  , and 

o o 
ID  ID 

< O 

II 

Equation  (4.20)  defines  a straight  line  through  U0  with  a step  length 
of  Ao  and  approaches  the  minimum  of  n-dimensional  space  in  the  direction 
of  the  vector  N.  The  rate  of  convergence  is  thus 


A 

n 

Ao 


(4.21) 


The  solution  is  obtained  when 


is  less  than  a given  small  number. 
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Obviously,  a lot  of  numerical  procedures  have  been  developed 
and  used  to  solve  general  nonlinear  problems.  Because  no  one  method 
appears  to  be  far  superior  to  all  the  others,  the  selection  of  a 
particular  technique  rests  on  the  characteristics  of  the  problem. 

For  a stress  concentration  problem,  it  is  presumed  that  the 
plastification  starts  from  the  regions  of  high  stress  gradients  and  stays 
in  these  areas  un  to  a certain  load  level.  Because  the  plastic  zones  are 
not  widely  spread  and  represent  only  a fraction  of  the  whole  structure,  a 
mildly  nonlinear  load-displacement  relation  is  anticipated.  To  solve 
this  type  of  problem,  a significant  amount  of  computational  effort  can 
be  saved  if  the  same  stiffness  is  used  in  the  iterative  method  to  obtain 
the  solution  corresponding  to  an  applied  increment  of  load.  The  structural 
stiffness  is  updated  only  on  request  to  speed  up  the  rate  of  convergence. 
The  incremental-iterative  method  [55],  based  on  a modified  Newton-Raphson 
procedure  is  thus  adopted  in  this  study.  This  method  is  usually  presented 
in  an  intuitive,  graphic  concept  as  shown  in  Fig.  10. 

To  avoid  the  accumulation  of  round-off  errors,  the  residual 
nodal  forces  from  the  previous  load  step  are  added  to  the  next  load 
increment. 

4 . 4 Outline  of  Numerical  Procedures 

Since  the  regions  with  highly  concentrated  stresses  are  modeled 
with  layered  3-D  isoparametric  elements,  it  seems  wasteful  to  examine  the 


f state  of  plastification  of  all  the  other  types  of  elements  if  those  elements 
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The  following  steps  are  carried  out  for  each  applied  load 

increment. 

(1)  Find  the  incremental  nodal  displacements  corresponding 
to  the  applied  load  increment  by  solving 

[ K] { AU } = { AP } 
n n 

(2)  Convert  the  incremental  displacements  to  strain  increments 
{Ae}n  and  then  use  the  elastic  material  properties  to 
calculate  a temporary  stress  increment,  {Aa'}n,  at  each 
integration  point  of  the  3-D  elements. 

(3)  Add  {Ao')n  to  {a)n  ^ and  compute  the  temporary  effective 
stress  o'  . 

Three  different  situations  may  occur  at  this  stage 

(4A)  If  c?  < a o , yielding  has  not  yet  happened  and  the 

temporary  stresses  are  the  actual  stresses.  Go  to  step  10. 

(4B)  If  on  1 >_a0,  yielding  is  already  occurring,  thus  the 

plastic  deformation  that  results  from  the  load  increment 
has  to  be  evaluated.  Find  [W]  as  defined  in  Eq.  (4.11) 
and  go  to  step  5. 

(4C)  If  an  i < a0  but  > a0  (see  Fig.  11) 

Yielding  begins  during  the  increment  when  the  stress 
vector  tries  to  penetrate  the  yield  surface.  The 
initial  yield  is  now  exceeded.  Therefore,  find  the 
intermediate  stress  value  by  multiplying  [ W]  by  a linear 
interpolation  factor  B [56]  which  is  defined  as 
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o0  - a 


8 = 1 - a = 1 - 


n-1 


o'  an  i 
n - n- 1 


Then,  proceed  to  step  5. 


-P 


(5)  Evaluate  the  effective  plastic  strain  increment  Ac 

1 1 

p 

and  the  plastic  strain  increment  of  each  component  {Ac  } 
by  employing  Eq.  (4.11)  and  Eq.  (4.8),  respectively. 

(6)  Find  the  actual  stress  increment  {Ao>n  by  substituting 
(Ae^  into  Eq.  (4.9),  then  add  to  the  previous  state  of 


stress. 


t0>n  = {a,n-l  * {4oln 


(7)  Compute  the  effective  stress  increment  Aa^,  then 
update  the  accumulated  effective  stress. 


°n  = Vl  + 4an 


(8)  Evaluate  the  unbalanced  nodal  forces. 

(9)  Update  [D  ] and  use  it  to  generate  the  new  element 

'-r 

stiffness  if  a new  structural  stiffness  matrix  is  required. 

(10)  Repeat  steps  2 to  9 until  all  3-D  elements  have  been 
examined. 

(11)  If  convergence  is  attained,  add  the  residual  nodal  forces 
to  the  next  load  increment.  Otherwise,  use  the  current 
force  unbalance  of  step  8 for  another  iteration  starting 
from  step  1 . 
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4. 5 Updating  Structural  Stiffness  Matrix 


When  a slow  rate  of  convergence  indicates  the  necessity  of 
updating  the  structural  stiffness,  usually  new  element  stiffness  matrices 
corresponding  to  the  current  material  properties  are  generated  for  all 
elements,  then  reassembled.  For  a stress  concentration  problem,  it  is 
presumed  that  those  elements  away  from  the  regions  of  rapid  variations 
in  displacements  are  still  in  their  elastic  range  even  up  to  the  load 
causing  severe  distress  in  the  intersection  region  of  the  structure. 

Thus,  it  is  more  efficient  to  use  the  difference  between  the  old  and 
the  present  element  stiffness  matrices  when  creating  the  new  system 
rather  than  reevaluating  and  reassembling  all  the  element  matrices. 

For  elements  that  have  yielded,  the  constitutive  law  gives 


{Aa}  = [Dep]{Ae> 


Therefore,  the  new  element  stiffness  matrix  is 


[K] 


[B]T[D  ][B]  dV 
vol  ep 


(4.22) 


(4.23) 


Equation  (4.23)  is  evaluated  at  integration  points.  For  those  points 
that  have  no  plastic  deformation,  the  [Dgp]  is  replaced  by  elastic 
material  property  matrix  [D]  as  described  in  Chapter  3. 

The  difference  or  change  in  the  element  stiffness  matrix  is 
calculated  as 

[AK]  - [K]„ew  - [K]0,d  (4.24) 

The  positions  of  [AK]  in  the  total  structural  stiffness  matrix  is 
determined  so  that  [AK]  can  be  added  to  the  old  structural  stiffness 
matrix  in  accordance  with  their  contributions  to  the  nodes  of  the 


structure.  To  save  storage  space,  only  the  lower  half  of  the  3-D 
element  matrix  is  stored  in  the  system. 

4.6  Evaluation  of  Excess  Nodal  Forces  [57] 

An  iterative  approach  is  commonly  employed  when  solving 
nonlinear  problems.  When  the  state  of  equilibrium  is  not  achieved, 
residual  stresses  exist  and  have  to  be  converted  to  unbalanced  nodal 
forces.  To  do  this,  the  residual  stresses  at  any  point  are  first 
evaluated  as  (Fig.  12) 

{Aaex)=  (Aa)g  - {Aa}^  (4.25) 


where  {Aalg  is  calculated  temporary  stress, 

(Ao)c  is  actual  stress  including  nonlinearity 

The  residual  nodal  forces  are  determined  by  converting  the  excess  nodal 
stresses  to  a system  loads  at  the  nodes  which  does  the  same  work  as  the 
excess  stresses  would  do  during  a virtual  displacement 


URC)  = l 


vol 


[B] '{Aaex} 


dV 


(4.26) 


Substituting  Eq.  (4.25)  into  Eq.  4.26,  we  get 

{ AR  } = £[  [ B]T{ Aa } R dV  - [B]T{ Aa}  d V 

c Jvol  B •'vol  c 


(4.27) 


The  substitution  of  Eq. (4. 8)  into  Eq.  (4.27)  yields 

B 

(AR  } = l\  [B]T{Aa}R  dV  - if  [B]T[D]({dd  - {dcP}  ciV  (4.28) 
Jvol  B •'vol  A 


If  the  initial  stiffness  is  used,  the  applied  load  increment  {AP}  can  be 
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expressed  as 


T B t 

{AP}  = l [B]  {Ao}r  dV  = l [B]  [D]{de}  dV  (4.29) 

Jvnl  B JV01  JA 


Therefore 


{AR  } = ][f  fB[B]T[D]{dcP}  dV 
c Jvol  JA 


(4.30) 


From  Eq.  (4.30),  it  is  seen  that  only  those  elements  which  have 
yielded  produce  excess  nodal  forces  and  need  to  be  included  when  integrating 
to  find  the  residual  forces.  This  saves  a lot  of  computer  operations  when 
solving  stress  concentration  problems  in  which  the  plastified  area  is  small 
compared  with  the  whole  structure. 

If  the  structural  stiffness  matrix  is  updated  (Fig.  13)  at 
point  E,  [D]  in  Eq.  (4.29)  is  replaced  by  a new  material  property  [Dgp] 
and  yields 


rB  T 

{AP}  = l [B]  [D-  ]{de}  dV 
Jvol  JA  ep 


where 


tV  = [D]  - [DP] 


fc  th  [Dp]  = CD]{|§}  [W] 


(4.31) 


[0]  - [Dep]  + [Dp] 


(4.32) 


i* 


Substituting  Eq.  (4.32)  into  (4.28),  we  get 


B B 

{AR  } = {AP}  - [ [B ]T [ D ] { de } dV  - xf  [B]T[Dl{dc}  dV 

c Jvol  M ep  Jvol  JA  p 


+ [B]T[D]{deP}  dV 

Jvol  M 


(4.33) 
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rB  t D fB  T 

(AR  } = l I [B]T[D]{dEP}  dV  - l [B]‘ [Dp]{dE}  dV  (4.34) 

c Jvol  JA  Jvol  JA 

The  first  term  in  Eq.  (4.34)  is  exactly  the  same  as  Eq.  (4.30). 
[Dp],  in  the  second  term,  is  zero  if  the  material  at  the  integration 
point  is  still  in  the  elastic  range  when  the  structural  stiffness  is 
updated.  Thus,  again,  only  the  plastified  elements  need  to  be  operated 
upon.  Obviously,  [Dp]  varies  from  point  to  point.  To  save  storage, 
only  {|^-}  is  kept  in  the  file  used  to  generate  [D  ]. 

dO  r 


57 


CHAPTER  5 

ELASTIC-PLASTIC  SOLUTION  OF  NORMALLY 
INTERSECTING  CYLINDERS 


5.1  General 

The  objective  of  the  present  study  is  to  develop  a finite 
element  method  to  solve  intersecting  cylinder  problems  including  in  that 
solution  any  elastic-plastic  behavior  that  might  develop.  Before  the 
selected  problem  was  investigated  several  numerical  examples  were  solved 
to  demonstrate  the  reliability  and  the  effectiveness  of  the  computational 
procedures.  Two  aspects  of  the  solution  process  had  to  be  evaluated. 

(1 ) Elastic  solution 

The  applicability  of  reduced  integration  techniques  to 
shell  structures  was  first  tested.  The  behavior  of  the 
elements  and  the  adequacy  of  the  proposed  discretization 
models  were  observed.  The  improvement  of  the  accuracy 
and  the  efficiency  of  the  procedures  suggested  the  use 
of  double  precision  and  the  secondary  storage  devices  in 
the  computer  work.  Before  a nonlinear  solution  was  sought 
an  elastic  solution  was  first  run.  The  regions  of  stress 
concentration  are  first  detected  from  these  elastic 
results.  Layered,  three-dimensional  elements  in  a finer 
mesh  are  then  applied  to  the  areas  with  high  stress 
gradients. 

(2)  Elastic-plastic  solutions 

The  incremental  loadings  are  applied  to  the  structure 
to  evaluate  the  corresponding  yielding  zones.  From  the 
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results,  the  progression  of  yield  zones  is  seen.  With 
the  availability  of  the  restart  feature  in  the  computer 
program,  the  structure  is  loaded  step  by  step  until 
failure. 

Finally,  a cyl inder-to-cyl inder  shell  intersection  was  modeled, 
subjected  to  internal  pressure,  and  solved  by  the  computer  program.  The 
elastic  solution  is  compared  with  experimental  data  and  Bakhrebah's 
finite  element  results.  The  plastic  solution  is  then  presented  and 
discussed. 

5.2  Elastic  Solutions 

5.2.1  Hyperbolic  Paraboloid  Shell 

The  applicability  of  the  reduced  integration  technique  to  a 
general  shell  element  is  demonstrated  by  solving  a hyperbolic  paraboloid 
shell  with  clamped  edges.  This  structure  is  subjected  to  c uniform  normal 
load.  The  geometrical  shape  of  the  midsurface  of  the  shell  is  described  by 

z = (f/ab)  xy 

where  f is  the  rise,  and  a and  b are  the  shell  spans  as  defined  in  Fig.  14. 

The  shell  is  antisymmetric  about  the  X and  the  Y axes. 

Therefore,  only  one  quarter  need  to  be  considered.  A 4 x 4 grid  was 
used  with  four  different  combinations  of  mesh  types  and  orders  of 
integration  rules. 

(1)  Regular  mesh  with  3x3  integration  points 

(2)  Regular  mesh  with  2x2  integration  points 
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(3)  Irregular  mesh  with  3x3  integration  points 

(4)  Irregular  mesh  with  2x2  integration  points 

The  vertical  displacement,  W,  across  midspan  (Y  = 0),  and  the  stress 

resultants  N , M at  integration  points  near  the  X axis,  are  shown  in 

y 

Figs.  15,  16,  and  17,  respectively.  From  a comparison  with  the  results 
obtained  by  Pecknold  and  Schnobrich  [58],  Choi  and  Schnobrich  [48], 
and  others,  several  points  bear  mentioning: 

(1)  All  four  of  the  cases  investigated  converged  to  give 

near  identical  solutions  for  W and  N . 

xy 

(2)  The  deviations  of  for  the  four  cases  considered  get 
larger  near  the  clamped  edge.  Actually,  My  in  both 
third  order  integration  schemes  "oscillated"  about  the 
results  obtained  by  other  researchers  while  the  second 
order  integration  procedures  gave  smooth  solutions. 

(3)  The  third  order  integration  rule  provided  slightly 
higher  results  in  the  regions  close  to  the  center  line 
but  smaller  values  near  the  clamped  edge  in  comparison 
with  the  results  obtained  from  second  order  integration. 

For  the  second  integration  rule,  a nonmonotonic  convergence 
is  anticipated. 

For  the  problem  studied  it  was  found  that  the  second  order 
integration  scheme  produced  better  results  than  did  the  higher  order 
integration.  This  was  true  not  only  for  displacements  but  more 
importantly  for  stress  resultants.  This  is  very  important  in  plastic 
analysis  when  the  integration  points  are  used  as  the  control  points 
to  predict  the  onset  of  plasticity  in  the  elements. 
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5.2.2  Pinched  Cylinder 


The  effectiveness  of  Ahmad's  element  when  used  with  a reduced 
integration  scheme  is  further  observed  in  the  solution  to  a pinched 
cylinder  problem.  This  problem  has  a stress  concentration  due  to  the 
application  of  point  forces.  A solution  to  the  pinched  cylinder  has 
been  published  by  Timoshenko.  This  problem  has  become  a base  against 
which  many  researchers  have  compared  their  solutions.  The  pinched 
cylinder  with  various  thickness  to  diameter  ratios  has  been  investigated 
by  Cantin  [59]  using  an  element  of  24  degrees-of-freedom.  Ashwell  and 
Sabir  [60]  also  solved  the  problem  by  employing  the  rectangular 
cylindrical  shell  element  with  20  degrees-of-freedom.  The  dimensions  of 
the  example  used  in  this  study  (Fig.  18)  were  taken  from  the  report  by 
Lindberg,  et  al.  [61].  This  shell  was  also  investigated  by  Razzaque  [62]. 
The  exact  solution  shown  for  comparison  was  obtained  by  Lindberg  using 
numerical  integration  of  Timoshenko's  equations. 

Because  of  symmetry,  only  one  octant  was  analyzed  and  that 
with  4x2  grids,  4 elements  in  the  hoop  direction  and  2 elements  in  the 
longitudinal  direction.  The  two  ends  of  the  cylinder  are  supported  by 
a diaphragm,  i.e.,  u = w = 0.  Equal  and  opposite  concentrate  loads  are 
applied  on  opposite  ends  of  the  diameter  at  the  center  points  of  the 
cyl inder. 


In  comparison  with  the  results 


provided  by  other  researchers. 


it  is  seen  (Figs.  19-22)  that  adequate  agreement,  in  many  cases  superior 
to  some  of  the  other  reported  results,  was  obtained  by  using  a relatively 
small  number  of  elements.  The  effectiveness  of  the  element  and  the 
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efficiency  of  the  reduced  integration  technique  make  nonlinear  analysis 
of  stress  concentration  problems,  such  as  the  intersecting  cylinder 
case  feasible. 


5.3  Elastic-plastic  Solutions 


5.3.1  Simply  Supported  Beam  and  Cantilevered  Beam 


As  a demonstration  of  the  adequacy  of  the  proposed  finite 
element  models  in  plastic  analysis,  two  elementary  but  representative 
sample  structures  were  solved.  A simply  supported  beam  and  a cantilevered 
beam  have  been  tested  with  the  same  structural  modeling,  but  different 
boundary  conditions.  The  structures  are  subjected  to  increasing  uniform 
loads.  An  elastic-ideal ly  plastic  material  behavior  is  assumed  in  both 
cases. 

A two-layered  system  with  3-D  quadrilateral  elements  was  first 
employed  in  the  region  of  high  bending  stresses,  i.e.,  the  fixed  end  of 
the  propped  cantilevered  beam  and  the  midspan  of  the  simple  beam. 
Transitional  elements  and  shell  elements  were  used  in  the  remaining 
portions  of  the  structures.  A second  order  integration  rule  was  used 
for  all  elements.  The  results  are  compared  with  the  theoretical  solutions 
provided  by  Prager  and  Hodge  [63]  in  Figs.  23 a and  24?.  In  these  figures. 
Wo  is  the  center  or  the  tip  deflection.  W$  is  the  center  or  the  tip 
deflection  at  the  maximum  load  for  which  the  structure  remains  entirely 
elastic,  and  p represents  the  nondimensional  load  parameter,  given  by 


P 


yield 


where  P is  applied  load  and  P0  = 4ba 


L and  t are  defined  as 
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shown  in  Figs.  23  and  24.  The  progression  of  the  elastic-plastic 
boundary,  through  the  cross  section,  is  quite  shallow  in  both  cases. 

A refined  mesh  with  four  layers  in  the  region  of  high  bending  stress 
was  then  adopted.  The  results  are  also  shown  in  Figs.  23a  and  24a. 

It  was  found  that  good  agreement  with  Hodge's  work  was  achieved. 

The  progression  of  the  plastic  zones  for  the  two  beams  are  shown  in 
Fig.  23b  and  Fig.  24b,  respectively.  The  flat  or  shallow  elastic- 
plastic  boundary  propagation  does  not  exist  for  most  of  the  intersecting 
structures  of  interest  in  the  present  study.  This  fact  is  taken  into 
account  when  the  modeling  is  selected  for  the  intersecting  cylinders. 

5.3.2  The  Thick-walled  Pressure  Vessel 

To  explore  the  general  elastic-plastic  behavior  of  cylinders 
under  increasing  internal  pressure,  an  infinitely  long  thick  hollow 
cylinder  with  internal  and  external  radii  of  lurnt  and  2unit,  respec- 
tively, has  been  analyzed.  This  problem  was  originally  investigated  by 
Hodge  [64]  using  a finite  difference  solution  with  the  case  of  elastic- 
plastic  materials.  Later,  both  Gupta  (23)  and  Salem  (55)  used  finite 
element  methods  to  obtain  comparison  solutions. 

A quarter  of  a unit-length  of  the  pressure  vessel  was  meshed 
by  8 three-dimensional  quadrilateral  elements  placed  in  four  layers 
through  the  thickness.  Symmetric  boundary  conditions  were  applied  to 
all  the  faces  except  the  internal  and  the  external  surfaces  of  the 
cylinder.  Second  order  reduced  integration  was  again  used  to  predict 
the  nonlinearity  of  material  and  to  evaluate  stresses.  The  results 


again  showed  a very  good  agreement  with  those  obtained  by  the  other 
researchers  (Figs.  25,  26  and  27). 

In  comparison  with  the  elastic  solution  at  p = 12.5  psi, 
it  is  concluded  that  stress  relief  occurs  in  the  hoop  stress  and  in 
the  axial  stress,  while  the  normal  stress,  o , is  almost  unaffected  by 
the  progression  of  the  yield  zone.  This  phenomenon  was  also  indicated 
by  Prager  [63].  If  the  pressurized  vessel  is  unloaded  elastically  after 
the  cylinder  has  been  stressed  into  the  plastic  range,  the  elastic  hoop 
stress  subtracts  from  the  actual  plastic  hoop  stress  to  give  a 
compressive  residual  stress  at  the  inside  wall.  This  compressive  stress 
is  found  very  useful  in  alleviating  high  stress  concentration  and  in 
industrial  applications  of  increasing  the  fatigue  life  of  highly 
pressurized  pipe  components.  From  the  results  at  other  loading  levels, 
it  is  found  that  the  plastic  zone  expands  in  accordance  with  the  peak  of 
the  hoop  stress  and  penetrates  through  the  cylindrical  wall  as  the  load 
is  increased.  When  the  structure  is  pressurized  to  12.5  psi,  the 
plasticity  propagates  to  a radius  of  1.4  times  the  internal  radius  as 
shown  in  Fig.  26b. 

5 . 4 Normally  Intersecting  Cylinders 
5.4.1  Introduction 

The  set  of  intersecting  cylinders  used  in  the  experiments 
conducted  at  Oak  Ridge  National  Laboratory  by  Corum  [7]  was  selected 
as  the  best  problem  to  be  analyzed  here  because  of  the  availability 
of  the  test  data.  The  problem  has  also  been  investigated  by  Prince 
and  Rashid  [18]  and  by  Bakhrebah  and  Schnobrich,  each  in  the  elastic 


range. 


64 


The  geometric  parameters  and  the  material  properties  of  the 
problem  are 

Nozzle-to-cyl inder  diameter  ratio,  d/D  = 0.5 
Cylinder  thickness-diameter  ratio,  T/D  = 100 
Nozzle  thickness-diameter  ratio,  t/d  = 100 
Nominal  hoop  stress  ratio,  s/S  = 1 
Outside  radius  of  the  cylinder,  R = 5 in. 

Cylinder  thickness,  T = 0.1  in. 

Modulus  of  elasticity,  E = 30,000,000  psi 
Poisson's  ratio,  v = 0.3 

Because  a mild  steel  material  is  being  simulated,  the  initial  yield  stress 
is  assumed  to  be  34  ksi  and  a very  shallow  hardening  modulus,  EH  = 700  ksi 
is  used.  The  geometrical  shape  of  the  intersection  curve  can  be  defined 
through  Eq.  (2.1)  by  setting  a = 0°  to  get 

x = r cos  4> 
y = r sin  4> 
z = - r2  cos2  cj> 

Substitution  of  the  internal  or  external  radii  of  the  nozzle 
and  the  cylinder  into  the  above  equations  produces  the  internal  or 
external  intersection  curve.  The  connection  of  the  outer  and  the  inner 
intersection  curves  determines  the  intersection  surface  between  the 
nozzle  and  the  cylinder. 

The  structure  is  subjected  to  increasing  internal  pressure. 

From  the  results  of  a limit  analysis,  the  lower  bound  [25]  and  upper 
bound  [26]  solutions  for  the  pressure  loads  to  initiate  yielding  and 
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to  produce  a mechanism  are  about  0.1  ksi  and  0.37  ksi , respectively. 

It  must  be  noted  that  the  upper  bound  0.37  ksi  is  for  perfect-plastic 
materials  without  considering  the  possibility  of  a local  failure  even 
though  a very  large  plastic  strain  may  have  occurred  at  some  point 
before  the  computed  failure  of  the  structure.  In  practical  applications, 
this  load  level  may  be  unattainable  because  the  large  deformations  may 
cause  the  structure  to  fail  physically  before  the  bound  is  reached. 
Besides,  the  load-displacement  curve  becomes  quite  shallow  when  the 
applied  load  approaches  this  limit.  This  creates  numerical  problems 
with  displacement  results  of  poor  accuracy  in  the  nonlinear  analysis. 
Therefore,  the  plastic  deformation  or  the  effective  plastic  strain  at 
some  critical  point  is  used  as  token  to  control  the  increase  of  the 
internal  pressure  in  this  study.  The  internal  pressure  loads  acting  on 
the  end  caps  of  the  vessel  and  the  nozzle  are  converted  to  axial  forces 
and  applied  to  the  walls  of  vessels. 

5.4.2  Discretization  Model  of  the  Intersecting  Cylinders 

The  structure  is  symmetric  about  the  XZ  and  the  YZ  planes 
(Fig.  28).  Thus  only  one-quarter  of  the  structure  need  be  considered. 
According  to  the  parametric  study  in  Chapter  2,  it  is  presumed  that 
the  stress  concentration  occurs  within  the  range  of  2.45/RT  from  the 
junction.  The  nozzle  and  the  cylinder  are  thus  modeled  with  layered 
3-D  elements  within  that  distance  from  the  intersection.  To  represent 
the  geometry  of  the  curved  intersection  surface,  quadrilateral  isopara- 
metric elements,  without  midside  nodes  in  the  thickness  direction,  are 
employed.  The  size  of  elements  is  increased  with  the  distance  away 
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from  the  junction.  The  thickness  of  the  vessels  is  divided  into  two 
layers  as  a result  of  consideration  of  the  following: 

(1)  The  thicknesses  of  the  vessels  are  very  thin  compared 
with  the  diameters.  Thus,  to  keep  the  depth  of  the  3-D 
elements  in  a proper  proportion  with  the  other  two 
dimensions,  a close  grid  refinement  through  the  thickness 
of  the  cylinder  is  not  desirable. 

(2)  The  half  bandwidth  of  the  total  structural  stiffness 
matrix  increases  greatly  with  the  number  of  layers. 

The  computer  time  is  proportional  to  the  total  degrees 
of  freedom  and  the  square  of  the  bandwidth. 

On  the  boundary  of  this  region,  3-D  transitional  elements  are  used  to 
connect  the  layered  elements  with  a single  element.  Shell  transitional 
elements  are  then  used  to  bridge  between  the  3-D  element  and  the  shell 
elements  which  are  used  throughout  the  remainder  of  the  structure  to 
reproduce  the  membrane  behavior  of  the  shell  structures. 

End  caps  are  used  at  the  ends  of  the  pipes,  and  the  boundary 
conditions  of  the  structure  are: 

At  the  end  cap  or  diaphragm  of  the  cylinder,  X = Z = B = 0 

At  the  end  cap  or  diaphragm  of  the  nozzle,  X = Y = 3 = 0 

On  the  symmetry  plane  YZ,  X = B = 0 

On  the  symmetry  plane  XZ,  Y = a = 0 

The  mesh  employed  in  this  analysis  (Fig.  28)  contains  507  nodes 
and  110  elements.  The  structure  has  about  2000  degrees-of-freedom  while 
the  half  bandwidth  spans  about  180  degrees-of-freedom.  The  input  data 
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was  generated  by  a computer  program.  The  intersecting  curve  function  was 
used  to  generate  the  coordinates  of  the  3-D  elements. 

The  structure  was  loaded  in  a step-by-step  manner.  At  the  end 
of  each  load  step,  the  generated  data  was  stored  on  a tape  for  restart  at 
a later  time.  The  structural  stiffness  matrix  was  updated  and  resolved 
twice  during  the  entire  solution  process.  Each  element  took  about  11 
seconds  of  IBM  360/75  computer  time  to  generate.  In  the  nonlinear  range, 
each  iteration  took  about  1 minute.  The  structure  was  loaded  by  eight 
loading  increments. 

The  reduced  integration  technique  was  also  employed  during  the 
whole  process.  The  second  order  integration  rule  was  used  for  all  the 
elements  to  evaluate  their  stiffness,  the  stresses  and  to  predict  the 
nonlinearity  of  the  structure. 


5.4.3  Elastic  Solution 


i 


V 


•i 


J 


The  elastic  solution  was  obtained  at  the  50  psi  level  of 
internal  pressure  load.  The  results  are  presented  and  compared  with 
the  experimental  values  and  with  the  finite  element  solutions  provided 
by  Bakhrebah.  These  comparisons  are  shown  in  Figs.  29  through  40.  Good 
agreement  is  observed  in  general  even  though  minor  deviations  appeared 
in  the  axial  stress  of  the  nozzle  near  the  0°  line  and  the  270°  line 
(Fig.  28).  The  axial  stresses  of  the  inside  surface  in  this  study 
converges  to  the  limiting  value  for  that  face,  that  is,  to  the  internal 
pressure.  This  is  expected  in  a proper  analysis  because  of  the 
satisfaction  of  the  equilibrium  condition  with  the  given  internal 
pressure  on  the  inner  surface. 
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From  the  elastic  solution,  several  observations  can  be  made: 

(1)  The  maximum  stresses  occur  at  the  junction  between  the 
nozzle  and  the  cylinder.  The  high  stresses  dissipate 
over  a range  less  than  1 . 5*^RT  from  intersection,  that 
means  about  1 in.  in  the  cylinder  and  0.5  in.  in  the 
nozzle  for  this  problem.  This  suggests  the  approximation 
2.45i/RT  is  on  the  conservative  side.  Away  from  this  region, 
both  the  hoop  stress  and  the  axial  stress  approach  the 
membrane  solution  values  and  stay  at  a constant  relation 
with  internal  pressure  P as 

oQ  = y-  = 2.5  ksi 
a = = 1-25  ksi 

(2)  The  hoop  stresses  in  the  nzzle  and  the  cylinder  near  the 
0°  line  are  in  tension  across  the  thickness  while  the  axial 
stresses  occur  as  tension  on  the  outside  fiber  and  compres- 
sion on  the  inside  fiber,  the  result  of  bending  moments. 

This  makes  the  inner  walls  highly  distorted  and  provides 
the  great  potential  for  yielding  to  initiate  here  when  the 
internal  pressure  is  increased. 

(3)  The  hoop  stress  on  the  cylinder  near  the  270°  line  results 
in  hoop  bending  to  balance  the  moment  produced  by  the 
internal  pressure.  Since  the  hoop  stress  in  this  region 

is  relatively  small  and  has  the  same  sign  as  the  longitudinal 
stress  both  on  the  inner  and  outer  fiber,  the  material  is 
thus  less  deformed  and  then  highly  resists  the  occurrence 

I; 

1 


69 


\ 


* 

* % 


f 

4* 


* 


of  plastification. 

5.4.4  Plastic  Solution 

The  initial  yield  load  obtained  from  this  study  is  about  0.11  ksi. 
This  value  compares  closely  with  the  result  of  0102  ksi  from  the  limit 
analysis  of  Ellyin  [25].  The  load-displacement  relations  for  points  A,  B 
and  C are  presented  in  Fig.  41.  The  stress  distributions  on  the  0°  line  of 
the  structure  at  load  level  P = 0.18  ksi  are  presented  in  Figs.  42  through 
49.  The  plastic  stress  distributions  are  compared  with  stress  values  that 
an  elastic  analysis  would  compute  at  the  same  load  level.  This  helps  in 
understanding  the  stress  development  as  the  structure  yields  and  visualizing 
the  residual  stresses  that  develop  upon  unloading. 

Figures  51  through  54  show  the  progression  of  the  plastified 
regions  in  the  structure  at  the  sections  containing  the  integration  points. 
The  plastic  action  occurs  as  the  internal  pressure  load  increases  from  0.1 
ksi  to  0.2  ksi.  The  boundaries  of  the  plastified  regions  were  drawn  by 
enveloping  the  plastified  integration  points  in  the  elements  at  the  various 
different  load  levels.  The  progression  of  the  yield  zones  along  the  inner 
and  the  outer  surfaces  of  the  structure  in  the  vicinity  of  the  intersection 
area  are  shown  in  Figs.  55  and  56,  respectively. 

From  the  numerical  results  obtained,  some  conclusions  have  been 

reached: 

(1)  Stress  relief  occurs  in  the  hoop  stress  in  the  plastified 
regions  while  the  axial  stresses  are  almost  unaffected  by 
the  progression  of  plastic  action.  The  stresses  converge 
to  the  membrane  solution  and  have  a constant  distribution 
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across  the  thickness  of  the  vessels  when  the  distance 
away  from  the  intersection  increases. 

(2)  The  initial  yield  occurs  at  the  inside  surface  of  the 
nozzle,  where  the  two  cylindrical  walls  intersect 
perpendicularly.  This  implies  that  the  bending  action 
resulting  from  the  discontinuity  of  this  region  requires 
special  consideration  in  the  structural  design. 

(3)  In  the  regions  around  the  intersection,  the  figures  show 
that  at  each  section  the  plastification  always  starts 
from  the  inside  layers,  then  gradually  penetrates  the 
walls  of  the  cylinders  from  the  inside  spreading  out  to 
the  outside. 

(4)  As  the  internal  pressure  is  increased,  the  plastified 
regions  progress  along  the  intersection  curve  and  spread 
out  following  the  distribution  of  high  stresses. 

(5)  The  elastic-plastic  boundaries  in  the  regions  of  the 
junction  are  parallel  to  the  longitudinal  direction  and 
layered  through  the  thickness,  while  away  from  these 
regions  the  boundaries  enveloping  the  plastified 
integration  points  become  steep.  This  implies  then 

that  membrane  axial  and  hoop  stresses  control  the  further 
progression  of  plastic  action.  Outside  the  immediate 
vicinity  of  the  intersection  the  normal  stress,  or>  is 
relatively  small  because  of  the  small  thickness-to- 
diameter  ratio.  Thus  the  normal  stress  does  not  play 
a significant  role. 
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The  largest  deformation,  about  0.01  in.,  occurred  at  point  C 
in  the  X direction  while  the  axial  displacement  of  the  main  vessel  (point  B) 
is  comparatively  small  and  almost  unaffected  by  the  plastification  of  the 
intersection.  The  effective  plastic  strain  at  the  inner  surface  point 
which  first  goes  plastic  increases  about  3 times  larger  than  initial  yield 
strain  when  the  pressure  load  is  0.2  ksi . Since  the  remaining  portions  out 
of  the  intersection  are  still  in  the  elastic  range,  an  increase  in  the 
thickness,  or  a fillet  reinforcement  around  the  junction,  may  be  adequate 
to  prevent  the  development  of  large  plastic  strains. 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  Conclusions 

A nonlinear  finite  element  method  has  been  presented  for 
analyzing  intersecting  cylinders.  From  the  sample  problems  investigated 
with  various  element  models  in  this  study,  it  is  found  that  good  results 
are  achievable  through  the  use  of  the  reduced  integration  technique.  When 
performing  the  plastic  analysis,  the  reduced  integration  points  serve  as 
the  control  points  to  predict  the  nonlinearity  of  the  material  behavior. 
Definition  of  the  progression  of  plastic  action  through  the  region  in  the 
vicinity  of  the  intersection  is  obtainable  with  the  incorporation  of 
layered  three-dimensional  elements.  Such  elements  are  used  in  the  regions 
of  high  stress  gradients  while  transitional  elements  and  two-dimensional 
curved  shell  elements  are  used  in  the  remaining  regions  to  economize  the 
computational  efforts. 

The  adequacy  of  the  discretized  model  of  the  intersecting 
cylinders  is  confirmed  by  comparing  the  elastic  solution  of  a normal 
cyl inder-to-cyl inder  intersection  with  experimental  data.  From  the 
nonlinear  solution  of  this  structure  it  is  found  that  yield  gradually 
penetrates  the  walls  of  the  cylinders  from  the  inside  face  to  the  outside 
face.  This  plastification  starts  from  the  region  where  the  two  cylindrical 
walls  intersect  perpendicularly,  then  spreads  out  along  the  intersection 
line  when  the  load  is  increased.  Around  the  intersection  line,  bending 
is  the  major  cause  of  yielding  while  away  from  this  region  membrane 
stresses  govern  the  progression  of  yielding. 
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Although  the  selected  intersecting  cylinders  fall  in  the 
category  of  thin  shell  structures,  the  procedures  can  be  applied  to 
general  cases  with  different  loading  patterns. 

6.2  Recommendations  for  Further  Studies 

6.2.1  Intersecting  Cylinders 

From  the  numerical  results  of  the  normally  intersecting  cylinders 
solved  in  this  study,  it  is  observed  that  the  plastified  region  is  limited 
to  a small  area  around  the  intersection  up  to  when  another  failure  mechanism 
forms  in  the  structure.  To  prevent  the  occurrence  of  large  plastic  strains 
in  this  area  or  to  keep  the  structure  in  a low  stress  field,  the  use  of 
reinforcement  or  a fillet  at  the  junction  is  of  particular  interest  in 
practical  design.  The  discretization  method  used  in  this  study  can  be 
utilized  to  investigate  the  effect  of  reinforcement  and  the  amount  it 
increases  the  strength  of  the  structure.  Research  can  also  be  expanded 
to  include  the  external  load  and  the  moment  caused  by  axial  loads  at  the 
end  of  the  branch  pipe. 

The  "K"  joint  is  another  type  of  connection  which  may  occur  in 
networks  of  pipes  or  offshore  oil-drilling  towers,  etc.  The  complicated 
geometrical  shape  of  intersection  curve  needs  special  considerations  in 
modeling  the  structure.  The  work  done  by  Greste  [18]  may  give  some 
general  ideas  for  the  start. 

6.2.2  Fatigue  Failure  of  Shell  Intersections 

Intersecting  cylinders  may  be  subjected  to  different  loading 
patterns  in  realistic  engineering  applications.  For  instance,  when  the 
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cyl inder-to-cyl inder  intersection  is  designed  to  operate  in  an  oceanic 
environment  the  loading  may  be  cyclic.  To  account  for  fatigue  behavior 
of  the  structure  becomes  important  in  this  case.  Since  elastic  unloading 
is  conducted,  the  Bauschinger  effect,  or  kinematic  hardening,  of  the 
materials  needs  to  be  considered.  Simplification  of  the  material  properti 
such  as  linear  strain  hardening  may  not  be  capable  of  representing  the 
realistic  situation  after  several  loading  cycles.  A more  specific 
description  of  the  material  properties  may  be  necessary.  The  work  done 
in  Ref.  65  may  be  helpful  in  understanding  the  general  behavior  in  this 
field. 

6.2.3  Fracture  Mechanics  of  Shell  Structures 

When  flaws  or  notches  exist  in  the  structure  under  certain 
unexpected  conditions,  fracture  of  the  material  may  occur  suddenly  causing 
the  leaking  or  the  complete  failure  of  the  structure.  Especially  when  the 
intersecting  cylinders,  with  their  high  stress  fields  around  the  inter- 
section, are  operated  in  an  environment  with  a temperature  lower  than  NDT 
of  the  material,  brittle  fracture  may  occur  without  any  warning  and  cause 
extensive  damage.  To  prevent  this  from  occurring  it  is  necessary  to 
investigate  the  fracture  mechanism  of  the  intersecting  cylinders  and  the 
progression  of  cracking.  A lot  of  attention  has  been  given  to  this  field. 
But,  the  application  of  the  finite  element  method  to  fracture  mechanics  is 
still  mainly  in  its  early  development  with  applications  in  two-dimensional 
plate  problems.  Very  little  has  been  publisned  on  the  application  of 
fracture  criteria  to  shell  structures.  There  are  two  major  difficulties: 

1.  The  lack  of  adequate  finite  element  models  to  represent 


75 


the  fracture  Dehavior  of  shell  structures. 

2.  A very  fine  mesh  around  the  tip  of  the  flaw  is  inevitable 
because  of  the  high  stress  concentrations  that  occur  in 
that  region. 

However,  the  layered  system  used  in  this  study  to  model  the  structure 
from  a fine  to  a coarse  mesh  around  the  tip  can  be  utilized  to  cut  down 
the  computational  expense. 
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Fig.  29  Cylinder-Cylinder  Intersection  Hoop  Stress  in 
the  Outside  Surface  of  Cylinder  near  0°  Line 


Fig.  30  Cylinder-Cylinder  Intersection  Hoop  Stress  in 
the  Inside  Surface  of  Cylinder  near  0°  Line 
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Fig.  31  Cylinder-Cylinder  Intersection  Axial  Stress  in 
the  Outside  Surface  of  Cylinder  near  0°  Line 
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Fig.  32  Cylinder-Cylinder  Intersection  Axial  Stress  in 
the  Inside  Surface  of  Cylinder  near  0°  Line 
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Fig.  42  Cylinder-Cylinder  Intersection  Hoop  Stress  in 
the  Outside  Surface  of  Cylinder  near  0°  Line 
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Fig.  43  Cylinder-Cylinder  Intersection  Hoop  Stress  in 
the  Inside  Surface  of  Cylinder  near  0°  Line 


Fig.  48  Cylinder-Cylinder  Intersection  Axial  Stress  in 
the  Outside  Surface  of  Nozzle  near  0°  Line 
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Fig.  49  Cylinder-Cylinder  Intersection  Axial  Stress  in 
the  Inside  Surface  of  Nozzle  near  0°  Line 
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APPENDIX  A 

THREE-DIMENSIONAL  ELEMENT 


A.  1 Elasticity  Matrix 

If  the  material  is  homogeneous  and  isotropic,  the  stress  components 
are  related  to  the  strain  components  by  an  equation  of  the  general  form. 

(a)  = [D]({e>  - (e0})  + {a0} 


where  {e0}  and  {o0}  are  the  initial  strain  and  initial  stress,  respectively. 
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A. 2 Stiffness  Matrix 


The  nodal  stiffness  matrix  [K.^]  can  be  expressed  as: 
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A. 3 T3D  Transformation  Matrix  [C] 


For  the  mid-node  i,  the  nodal  displacement  can  be  transformed  to  top 
and  bottom  nodes  by 
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APPENDIX  B 

AHMAD'S  SHELL  ELEMENT 


B.l  Shape  Function 


Ni 

0 

0 

N. 

v 

2 

en 

-Ni 

~r 

012 

[N*]i  = 

0 

Ni 

0 

Ni 

V 

2 

021 

-Ni 

V 

2 

022 

0 

0 

Ni 

Ni 

V 

~r 

031 

•Ni 

V 

2 

032 

where  [6]  is  defined  as 
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B.2  Material  Properties  [D * ] 

The  [D]  matrix  for  an  anisotropic  material  is  constructed  i 
coordinate  system  as: 
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where  k is  taken  as  1.2,  as  indicated  by  Zienkiewicz,  to  construct  the 
same  shear  strain  energies  as  in  real  distributions  through  the  thickness 
of  shells. 


B.3  [B1]  Matrix  for  Ahmad's  Element 
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where 


[9]  = [vp  v2,  v3] 
[J*]  = [J]"1  , and 
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The  [B1]^  for  node  i can  be  expressed  as: 
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B.4  Shell  Transition  Element 


When  the  departure  concept  is  incorporated  in  the  two-dimensional 
curved  shell  element,  the  shape  functions  of  the  nodes  that  are  intended 
to  be  converted  then  connected  to  a 3D  element  can  be  written  as: 

for  corner  nodes,  = -1,  = ±1 
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for  midside  nodes,  - -1,  = 0 
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The  [B1]  matrix  can  be  obtained  by  substituting  these  new  shape 
functions  into  Eq.  (3.29),  then  use  [B * ] to  construct  a modified  stiffness 
matrix  [Kj  1 . The  stiffness  matrix  of  the  shell  transition  element  is 
obtained  by 


[kt]  = [cs]'[k;][cs] 


where  [C$]i  for  node  i is  defined  as 
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APPENDIX  C 


VON  MISES  CRITERION 


From  the  experimental  investigations  [50],  it  has  been  demonstrated 
that  hydrostatic  stress  does  not  produce  yielding  in  a ductile  material 
even  for  very  high  hydrostatic  stress.  Based  on  this  fact,  the  distortion 
energy  theory,  or  the  Von  Mises  yield  criterion,  assumes  that  yielding 
begins  when  the  distortion  energy  equals  the  distortion  energy  at  yield 
in  the  simple  tension.  The  distortion  energy  of  a material  can  be 
written  in  general  form  as: 


ii  - — 1 - 3 2 

ud  " 2G  4G  cct 


(C.l) 


2[TkT 


For  simple  tension,  the  distortion  energy  is 


Ud  2G  = 2G  °° 


(C  .2) 


a o = first  initial  yield  stress  of  the 
simple  tension 


From  Eq.  (C.l)  and  Eq.  (C.2),  we  get 


J2  = T Oo 


(C.  3) 


Since  J = j ToCt’  there'fore  Eq.  (C.3)  can  be  written  as 
Toct  “ Toct,  simple  tension 


(C.4) 


Equation  (C.4)  indicates  that  yielding  occurs  when  the  octahedral  shear 
stress  reaches  the  octahedral  shear  stress  at  yield  in  simple  tension. 
The  Von  Mises  yield  criterion  usually  provides  good  agreement  with 
experimental  data  if  the  material  is  homogeneous  and  isotropic. 
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